使用 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.