使用 Matlab 的 bvp4c 求解边值问题

Kuiken 利用相似变换，得到如下非线性微分方程

满足如下边界条件

其中， 表示对 求导，为普朗特数. 此方程是耦合的非线性边值问题，在无穷远点具有奇性.

 

当 时，使用Matlab的bvp4c求解如下：

将原方程转化为一阶方程组

% kuikenode.m
function df=kuikenode(eta,f)
sigma=1; 
df=[ f(2)
     f(3)
     f(2)^2-f(4)
     f(5)
     3*sigma*f(2)*f(4)];

输入边界条件

% kuikenbc.m
function res=kuikenbc(f0,finf)
res =[f0(1)
      f0(2)
      f0(4)-1
      finf(2)
      finf(4)];

以常数作为初始猜测解

% kuikeninit.m
function v=kuikeninit(eta)
v =[ 0
     0
     1
     0
     0];

调用bvp4c求解，注意此处无界区间被截断

% solve.m
clc;
clear;
infinity=30;
solinit=bvpinit(linspace(0,infinity,5),@kuikeninit);
options=bvpset('stats','on','RelTol', 1e-12);
sol=bvp4c(@kuikenode,@kuikenbc,solinit,options);
eta=sol.x;
g=sol.y;
fprintf('\n');
fprintf('Kuiken reports f''''(0) = 0.693212.\n')
fprintf('Value computed here is f''''(0) = %7.7f.\n',g(3,1))
fprintf('Kuiken reports %c''(0) = -0.769861.\n', char([952]))
fprintf('Value computed here is %c''(0) = %7.7f.\n',char([952]),g(5,1))

clf reset
subplot(1,2,1);
plot(eta,g(2,:));
axis([0 infinity 0 1]);
title('Kuiken equation, \sigma =1.')
xlabel('\eta')
ylabel('df/d\eta')

subplot(1,2,2);
plot(eta,g(4,:));
axis([0 infinity 0 1]);
title('Kuiken equation, \sigma = 1.')
xlabel('\eta')
ylabel('\theta')
shg

运行如下：

 

 

The solution was obtained on a mesh of 105 points.
The maximum residual is 9.866e-013. 
There were 6332 calls to the ODE function. 
There were 260 calls to the BC function. 


Kuiken reports f''(0) = 0.693212.
Value computed here is f''(0) = 0.6932116.
Kuiken reports θ'(0) = -0.769861.
Value computed here is θ'(0) = -0.7698611.