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