常微分方程组边值问题解法打靶法Shooti ng Method (shoot in g.m )% 丁靶法求常微分方程的边值问题function [x,a,b ,n]=shooti ng(fu n, xO,x n, eps) if nargin<3eps=1e-3;endx1=x0+ra nd;[a,b]=ode45(fu n, [0,10],[0,x0]');c0=b(le ngth(b),1);[a,b]=ode45(fu n, [0,10],[0,x1]');c1=b(le ngth(b),1);x2=x1-(c1-x n)*(x1-x0)/(c1-c0);n=1;while (no rm(c1-x n)>=eps & no rm(x2-x1)>=eps) x0=x1;x 仁x2;[a,b]=ode45(fu n,[ 0,10],[0,x0]');cO=b(le ngth(b),1);[a,b]=ode45(fu n,[ 0,10],[0,x1]');c1= b(le ngth(b),1)x2=x1-(c1-x n)*(x1-x0)/(c1-c0);n=n+1;endx=x2;应用打靶法求解下列边值问题:y 10 0 解：将其转化为常微分方程组的初值问题命令:xO=[O:O.1:1O];y0=32*((cos(5)-1)/si n( 5)*si n(x0/2)-cos(x0/2)+1); plot(xO,yO,'r')hold on [x,y]=ode45('odebvp',[0,10],[0,2]'); plot(x,y(:,1))[x,y]=ode45('odebvp',[0,10],[0,5]'); plot(x,y(:,1))[x,y]=ode45('odebvp',[0,10],[0,8]'); plot(x,y(:,1)) [x,y]=ode45('odebvp',[0,10],[0,10]'); plot(x,y(:,1))dy i dx y 2 dy 2 dx y i 0y 4 y odydx X0真实解30'12^4567^9 10函数：(odebvp.m)泌值常微分方程(组)函数function f=odebvp(x,y)f(1)=y(2)；f(2)=8-y(1)/4;f=[f(1)；f(2)]；命令：[t,x,y, n]=shoot in g('odebvp',10,0,1e-3)计算结果：(eps=0.001 )t=11.9524plot(x,y(:,1))x0=[0:1:10];y0=32*((cos(5)-1)/si n( 5)*si n(x0/2)-cos(x0/2)+1); hold onplot(xO,yO, ' o')有限差分法 Finite Differenee Methods FDM同上例:Y i i2Y i y i ih 2 2 ——14函数：(differe nce.m )%有限差分法求常微分方程的边值问题 function [x,y]=differe nce(xO,x n,yO,yn,n) h=(x n-xO)/n;a=eye( n-1)*(-(2-h A 2/4)); for i=1: n-2 a(i,i+1)=1; a(i+1,i)=1; endb=o nes( n-1,1)*8*hA2; b(1)=b(1)-0; b(n-1)=b( n-1)-0; yy=a\b;x(1)=x0;y(1)=y0; for i=2: n x(i)=x0+(i-1)*h;Y i i 2若划分为10个区间，则：h- Y i Y i 1 8h 24(differe nce.m )1h 2 2 411Y 1 Y 28h 28h 20 h 2Y n 2 8h 224 1Y n 18h 2 0.2h 124y(i)=yy(i-i);endx(n)=xn;y(n)=yn;命令：[x,y]=differe nce(0,10,0,0,100);计算结果：xO=[O:O.1:1O];y0=32*((cos(5)-1)/si n( 5)*si n(x0/2)-cos(x0/2)+1);plot(xO,yO,'r')hold on[x,y]=differe nce(0,10,0,0,5);plot(x,y,'.')[x,y]=differe nce(0,10,0,0,10);plot(x,y,'--')[x,y]=differe nce(0,10,0,0,50);plot(x,y,'-.')正交配置法Orthogonal Collocatioin Methods CM构造正交矩阵函数(collmatrix.m )%E交配置矩阵(均用矩阵法求对称性与非对称性正交配置矩阵) function [am,bm,wm,an,bn,wn]=collmatrix(a,m,fm,n,fn) 真实解的3 4 6 7 9 9 1DXxO=symm(a,m,fm); %a为形状因子;m为零点数;fm为对称的权函数(0为权函数1,非0为权函数1-x A2)for i=1:mxm(i)=x0(m+1-i);endxm(m+1)=1;for j=1:m+1for i=1:m+1qm(j,i)=xm(j)A(2*i-2);cm(j,i)=(2*i-2)*xm(j)A(2*i-3);dm(j,i)=(2*i-2)*(2*i-3+(a-1))*xm(j)A(2*i-3+(a-1)-1-(a-1));endfmm(j)=1/(2*j-2+a);endam=cm*i nv(qm);bm=dm*i nv(qm);wm=fmm*i nv(qm);x1=unsymm(n,fn); %n 为零点数;fn为非对称的权函数(0为权函数1,非0为权函数1-x) xn (1)=0; for i=2: n+1xn( i)=x1( n+2-i);endxn(n+2)=1;for j=1: n+2for i=1: n+2qn (j,i)=x n(j)A(i-1);if j==0 | i==1cn (j,i)=0;elsecn (j,i)=(i-1)*x n(j)A(i-2);endif j==0 | i==1 | i==2dn (j,i)=O;elsedn (j,i)=(i-2)*(i-1)*x n(jF(i-3);endendfnn (j)=1/j;endan=cn *i nv(qn);bn=dn *i nv(qn); wn=fnn*inv(qn); %E交多项式求根(适用于对称问题)fun ctio n p=symm(a,m,fm) %a for i=1:mc1=1;c2=1;c3=1;c4=1;for j=0:i-1c1= c1*(-m+j);if fm==0c2=c2*(m+a/2+j);% elsec2=c2*(m+a/2+j+1);% endc3=c3*(a/2+j);c4=c4*(1+j);endp(m+1-i)=c1*c2/c4/c3;endp(m+1)=1;%为多项式系数向量p=sqrt(roots(p)); 为形状因子,m为配置点数,fm为权函数权函数为1权函数为1-x A2,求出根后对对称问题还应开方才是零点%E交多项式求根（适用于非对称性问题）解：（1）标准化令x r/R ，y C/C S 代入微分方程及边界条件得：丄2 x 2dy36yx 2dx dxx 0単0 dx x 1, y 1fun ctio n p=un sym m(n,fn) if fn==0r(1)=(-1)A n*n *( n+1);% 权函数为 1elser(1)=(-1)A n*n*(n+2);% 权函数为 1-x end for i=1: n-1 if fn==0r(i+1)=( n-i)*(i+n+1)*r(i)/(i+1)/(i+1);% elser(i+1)=( n-i)*(i+n+2)*r(i)/(i+1)/(i+1);% end end for j=1: np( n+1-j)=(-1)A(j+1)*r(j); endp(n +1)=(-1)A( n+1); p=roots(p);权函数为1权函数为1-x 应用正交配置法求解以下等温球形催化剂颗粒内反应物浓度分布，其浓度分布的数学 模型为：丄r 2 drr r2dC r dr36C R 2dC 0，丁 dr(2)离散化N 1B ji y i 36y j 0i 1j 1, 2,N 1(3)转化为代数方程组(以 N 3为例)因为y N i 目41，所以整理上式得:B 11 36 B 12 B 13B 21 B 22 36 B 23 B 31 B 32 B 33 36B 41 B 42 B 43y iy 2 y 3 B 14 B 24B 34 B 44 36 本例中的代数方程组为线性方程组，可采用线性方程组的求解方法；若为非线性方程组则 采用相应的方法求解。

命令：N=3,权函数为1-x [am,bm,wm,a n,bn,wn ]=collmatrix(3,3,1,3,1);b1=bm;for i=1:4 b1(i,i)=bm(i,i)-36; end a0=b1(1:4,1:3); b0=-b1(1:4,4); y=aO\bO; y(4)=1;p=exam31(3,3);(注意要对文件修改权函数为 x=[0.3631,0.6772,0.8998,1]; %零点 plot(x,y,'o')hold onx0=0:0.1:1; % 真实解 y0=si nh(6*x0)./x0/si nh(6); plot(xO,yO,'r') (只用对称性配置矩阵) 1-x 2)若权函数改为1,则以下语句修改，其他不变 [am,bm,wm,a n,b n,w n]=collmatrix(3,3,0,3,1);(只用对称性配置矩阵)r1 2 3 4 yyyy3^142434BBB33433421112 3 4BBBp=exam31(3,3);(注意要对文件修改权函数为1)x=[0.4058,0.7415,0.9491,1]; % 零点计算结果:权函数为1- x权函数为1%11正交配置法0.9真实解0.80.7 ・0.6 -y 0.5 -0.4 - -0.3 -X：0.2 --0.1 --1ii 1 ■0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1x边值问题的MatLab解法精确解:2x函数：(collfuni.m ) fun ctio n f=collfu ni(x,y) f(i)=y(2);f(2)=4*y(1);f=[f(1);f(2)];(collbci.m ) function f=collbc1(a,b) f=[a(1)-1;b(1)-exp(2)];命令：soli ni t=bvpi ni t([0:0.1:1],[1,1]) sol=bvp4c(@collfu n1,@collbc1,soli nit) plot(sol.x,sol.y)hold onplot(sol.x,exp(2*sol.x),'*')yi y2 真实解y 4y 0 x 1y 0 1, y 1 e2y i yy2 4y iy i 0 i, y i i e2真实解精确解:函数：(collfun2.m )function f=collfun2(x,y) f(1)=y(2)；f(2)=(1-x.A 2).*exp(-x)+2*y(1)-(x+1).*y(2); f=[f(1)；f(2)]；(collbc2.m )function f=collbc2(a,b) f=[a(2)-2;b(2)-exp(-1)];命令：soli nit=bvpi ni t([0:0.1:1],[1,1]);sol=bvp4c(@collfu n2,@collbc2,soli nit); plot(sol.x,sol.y) hold onplot(sol.x,(sol.x-1).*exp(-sol.x),'*')y x 1 y 2y 1 x 2 e x0 x 1y 0 2, y 1 1/ey i y 2y 2 1 x 2 e x 2y 1 % 0 2, y 1 1x 1 y 21/e 真实解函数：(collfu n3.m ) function f=collfun3(x,y) f(i)=y(2);f(2)=(2-log(x))./x+y(1)./x-y(2).A2; f=[f(i)；f(2)]；(collbc3.m )function f=collbc3(a,b) f=[2*a(i)-a (2);b(2)-1.5];命令：soli nit=bvpi ni t([1:0.1:2],[1,1]);sol=bvp4c(@collfu n3,@collbc3,soli nit); plot(sol.x,sol.y) hold onplot(sol.x,sol.x+log(sol.x),'*')真实解2y yy12ln x 1 x 2x x2y 1y 1 0,y 23/2精确解：y x In xy 2 2 In xx2y i 1 y 2 i 0,y i y 2xy 2 23/2在260 C 的基础面上，为促进传热在此表面上 增加纯铝的圆柱形肋片，其直径为25mm 高为150mm 该柱表面受到16 C 气流的冷却，气流与肋 片表面的对流传热系数为15W/m 2 K ，肋端绝热;肋片的导热系数为 236W/m K ，假设肋片的导热 热阻与肋片表面的对流传热热阻相比可以忽略；试 求肋片中的温度分布，及单个肋片的散热量为多 少？解：根据以上条件可知：导热热阻与对流热阻相比可以忽略，则在肋片径向上没有温 度分布，在轴向上存在温度分布，外界气流与肋片的对流传热则可转化为内热源；故该问 题为导热系数为常数的一维稳定热传导，其导热微分方程为：这是两点边值的常微分方程求解问题，故可转化为如下形式:y iy 2hp y yA Cy 2 H 0, 0边界条件为：x 0时，t 0 分析解：t t t 0 td 2t dx 2hp t t A C260 C （肋根）； x H 时，dtdx0 （肋端绝热）。