使用 Matlab 的 bvp4c 求解边值问题

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

满足如下边界条件

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

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

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

% kuikenode.mfunction 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.mfunction res=kuikenbc(f0,finf)res =[f0(1)      f0(2)      f0(4)-1      finf(2)      finf(4)];

以常数作为初始猜测解

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

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

% solve.mclc;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 resetsubplot(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.