用Matlab bvp4c 解带未知参数的常微分方程边值问题,怎么给出合适的初值?



用算出来的一个结果做为初值,即continuation 方法。下面用具体的例子作为说明。例子来自Solving ODEs with Matlab

方程为:

\begin{equation*} \begin{split} f'''-R[(f')^2-ff'']+RA=&0\\ R''+Rfh'+1=&0\\ \theta ''+P_ef\theta '=&0 \end{split} \end{equation*}

边界条件为:

\begin{equation*} \begin{split} f(0)=f'(0)=&0\\ f(1)=f'(1)=&1\\ h(0)=h(1)=&0\\ \theta(0)=&0\\ \theta(1)=&1 \end{split} \end{equation*}

方程描述流体描述流体流过竖直管道的问题,其中 $R$ 为常数,在模型中为雷诺数。常数$P_e=0.7R$,在模型中为佩克莱特数。$A$ 为我们要计算的未知常数。

给定 $R=100, 1000, 10000$,分别解方程。$R=100$时,很容易得到结果,$R=1000, 10000$时程序不收敛,可以把$R=100$时的结果作为初值代入,则顺利得到结果。

下面给出两个程序。

不用雅克比矩阵

程序见 https://ww2.mathworks.cn/matlabcentral/fileexchange/3819-tutorial-on-solving-bvps-with-bvp4c?focused=6778130&tab=function

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
function ex8bvp ()
%EX8BVP Example 8 of the BVP tutorial.
% These are equations describing fluid injection through one side of
% a long vertical channel, considered in Example 1.4 of U.M. Ascher,
% R.M.M. Mattheij, and R.D. Russell, Numerical Solution of Boundary
% Value Problems for Ordinary Differential Equations, SIAM, 1995
%
% The differential equations
%
% f''' - R*[(f')^2 - f*f''] + R*A = 0
% h'' + R*f*h' + 1 = 0
% theta'' + P*f*theta' = 0
%
% are to be solved subject to boundary conditions
%
% f(0) = f'(0) = 0
% f(1) = 1
% f'(1) = 0
% h(0) = h(1) = 0
% theta(0) = 0
% theta(1) = 1
%
% Here R and P are known constants, but A is determined by the boundary
% conditions.
%
% For a Reynolds number R = 100, this problem can be solved with crude
% guesses, but as R increases, it becomes much more difficult because of
% a boundary layer at x = 0. The solution is computed here for R = 10000
% by continuation, i.e., the solution for one value of R is used as guess
% for R = R*10. This is a comparatively expensive problem for BVP4C because
% a fine mesh is needed to resolve the boundary layer and there are 7
% unknown functions and 1 unknown parameter.
<p>% The solution is first sought for R = 100.<br />
R = 100;</p>
<p>options = []; % place holder</p>
<p>% A crude mesh of 10 equally spaced points and a constant 1<br />
% for each solution component is used for the first value of R.<br />
% The unknown parameter A is guessed to be 1.<br />
solinit = bvpinit(linspace(0,1,10),ones(7,1),1);</p>
<p>sol = bvp4c(@ex8ode,@ex8bc,solinit,options,R);</p>
<p>fprintf('For R = %5i, A = %4.2f.\n',R,sol.parameters);<br />
figure<br />
lines = {'k-.','r--','b-'};<br />
plot(sol.x,sol.y(2,:),lines{1});<br />
axis([-0.1 1.1 0 1.7]);<br />
title('Fluid injection problem');<br />
xlabel('x');<br />
ylabel('f'' (x)');<br />
drawnow</p>
<p>% The solution is computed for larger R by continuation, i.e.,<br />
% the solution for one value of R is used as guess for the next.<br />
hold on<br />
for i=2:3<br />
R = R*10;<br />
sol = bvp4c(@ex8ode,@ex8bc,sol,options,R);<br />
fprintf('For R = %5i, A = %4.2f.\n',R,sol.parameters);<br />
plot(sol.x,sol.y(2,:),lines{i});<br />
drawnow<br />
end<br />
legend('R = 100','R = 1000','R = 10000',1);<br />
hold off</p>
<p>% --------------------------------------------------------------------------</p>
<p>function dydx = ex8ode(x,y,A,R);<br />
%EX8ODE ODE function for Example 8 of the BVP tutorial.<br />
P = 0.7<em>R;<br />
dydx = [ y(2)<br />
y(3)<br />
R</em>(y(2)^2 - y(1)<em>y(3) - A)<br />
y(5)<br />
-R</em>y(1)<em>y(5) - 1<br />
y(7)<br />
-P</em>y(1)*y(7) ];</p>
<p>% --------------------------------------------------------------------------</p>
<p>function res = ex8bc(ya,yb,A,R)<br />
%EX8BC Boundary conditions for Example 8 of the BVP tutorial.<br />
res = [ya(1)<br />
ya(2)<br />
yb(1) - 1<br />
yb(2)<br />
ya(4)<br />
yb(4)<br />
ya(6)<br />
yb(6) - 1];</p>

显示给出雅克比矩阵

程序见Solving ODEs with Matlab 一书 pp192

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
function ch3ex6
% CH3EX6 example 3.5.6 of book Solving ODEs with Matlab
<p>global R<br />
color = [ 'k', 'r', 'b' ];<br />
options = bvpset('FJacobian',@Jac,'BCJacobian',@BCJac,...<br />
'Vectorized','on');<br />
R = 100;<br />
sol = bvpinit(linspace(0,1,10),ones(7,1),1);<br />
hold on<br />
for i = 1:3<br />
sol = bvp4c(@ode,@bc,sol,options);<br />
fprintf('For R = %5i, A = %4.2f.\n',R,sol.parameters);<br />
plot(sol.x,sol.y(2,:),color(i));<br />
axis([-0.1 1.1 0 1.7]);<br />
drawnow<br />
R = 10<em>R;<br />
end<br />
legend('R = 100','R = 1000','R = 10000');<br />
hold off<br />
%================================================================<br />
function dydx = ode(x,y,A);<br />
global R<br />
P = 0.7</em>R;<br />
dydx = [ y(2,:); y(3,:); R<em>(y(2,:).^2- y(1,:).</em>y(3,:)-A);...<br />
y(5,:); -R<em>y(1,:).</em>y(5,:)-1; y(7,:); -P<em>y(1,:).</em>y(7,:) ];<br />
function [dFdy,dFdA] = Jac(x,y,A)<br />
global R</p>
<p>dFdy = [ 0, 1, 0, 0, 0, 0, 0<br />
0, 0, 1, 0, 0, 0, 0<br />
-R<em>y(3), 2</em>R<em>y(2), -R</em>y(1), 0, 0, 0, 0<br />
0, 0, 0, 0, 1, 0, 0<br />
-R<em>y(5), 0, 0, 0, -R</em>y(1), 0, 0<br />
0, 0, 0, 0, 0, 0, 1<br />
-7/10<em>R</em>y(7), 0, 0, 0, 0, 0, -7/10<em>R</em>y(1) ];<br />
dFdA=[0;0;-R;0;0;0;0];<br />
function res = bc(ya,yb,A)<br />
res = [ ya(1); ya(2); yb(1)-1; yb(2);...<br />
ya(4); yb(4); ya(6); yb(6)-1 ];<br />
function [dBCdya,dBCdyb,dBCdA] = BCJac(ya,yb,A)<br />
dBCdya=[1,0,0,0,0,0,0<br />
0, 1, 0, 0, 0, 0, 0<br />
0, 0, 0, 0, 0, 0, 0<br />
0, 0, 0, 0, 0, 0, 0<br />
0, 0, 0, 1, 0, 0, 0<br />
0, 0, 0, 0, 0, 0, 0<br />
0, 0, 0, 0, 0, 1, 0<br />
0, 0, 0, 0, 0, 0, 0 ];<br />
dBCdyb=[0,0,0,0,0,0,0<br />
0, 0, 0, 0, 0, 0, 0<br />
1, 0, 0, 0, 0, 0, 0<br />
0, 1, 0, 0, 0, 0, 0<br />
0, 0, 0, 0, 0, 0, 0<br />
0, 0, 0, 1, 0, 0, 0<br />
0, 0, 0, 0, 0, 0, 0<br />
0, 0, 0, 0, 0, 1, 0 ];<br />
dBCdA = zeros(8,1);</p>

标签: matlab, bvp4c, 常微分方程, 边值问题

添加新评论

captcha
请输入验证码