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

\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*}

## 不用雅克比矩阵

 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>

## 显示给出雅克比矩阵

 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>