常微分方程组边值

1、常微分方程组边值问题解法打靶法 Shooting Method (shooting.m )%打靶法求常微分方程的边值问题function x,a,b,n=shooting(fun,x0,xn,eps)if nargin=eps & norm(x2-x1)=eps)x0=x1;x1=x2;a,b=ode45(fun,0,10,0,x0);c0=b(length(b),1);a,b=ode45(fun,0,10,0,x1);c1=b(length(b),1)x2=x1-(c1-xn)*(x1-x0)/(c1-c0);n=n+1;endx=x2;应用打靶法求解下列边值问题：,248工dx 4y 00

2、y 100解：将其转化为常微分方程组的初值问题y2y4 y。tdy1dx仪总8dx30dydx x0命令:x0=0:0.1:10;真实解y0=32*(cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1);plot(x0,y0,r)hold onx,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)函数：

3、(odebvp.m )%边值常微分方程(组)函数function f=odebvp(x,y)f(i)=y(2);f(2)=8-y(1)/4;f=f(1);f(2);命令：t,x,y,n=shooting(odebvp,10,0,1e-3)计算结果：(eps=0.001)t=11,9524plot(x,y(:,1)x0=0:1:10;y0=32*(cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1);hold onplot(x0,y0, o)FDM (difference.m )有限差分法 Finite Difference Methods同上例:.2d y 8 3yi

4、 i 2v yi i 8 至dx24h24c h22x 12yiyi 1 8h4若划分为10个区间，则:h242 h-2V18h2 0y28h2Yn 28h22Yn 18h 0h2函数：(difference.m )%有限差分法求常微分方程的边值问题function x,y=difference(x0,xn,y0,yn,n)h=(xn-x0)/n;a=eye(n-1)*(-(2-hA2/4);for i=1:n-2a(i,i+1)=1；a(i+1,i)=1;endb=ones(n-1,1)*8*hA2;b(1)=b(1)-0;b(n-1)=b(n-1)-0;yy=ab;x(1)=x0;y(1)

5、=y0;for i=2:nx(i)=x0+(i-1)*h;y(i)=yy(i-1);endx(n)=xn;y(n)=yn;命令：x,y=difference(0,10,0,0,100);计算结果：x0=0:0.1:10;y0=32*(cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1);真实解plot(x0,y0,r)hold onx,y=difference(0,10,0,0,5);plot(x,y,.)x,y=difference(0,10,0,0,10);plot(x,y,-)x,y=difference(0,10,0,0,50);正交配置法 Orthogona

6、l Collocatioin MethodsCM构造正交矩阵函数(collmatrix.m )%正交配置矩阵(均用矩阵法求对称性与非对称性正交配置矩阵)function am,bm,wm,an,bn,wn=collmatrix(a,m,fm,n,fn)x0=symm(a,m,fm); %a为形状因子;m为零点数;fm为对称的权函数(0为权函数1,非0为权函数1-xA2)for i=1:mxm(i)=x0(m+1-i);endxm(m+1)=1;for j=1:m+1for i=1:m+1 cm(j,i)=(2*i-2)*xm(j)A(2*i-3);dm(j,i)=(2*i-2)*(2*i-3+

7、(a-1)*xm(j)A(2*i-3+(a-1)-1-(a-1);endfmm(j)=1/(2*j-2+a);endam=cm*inv(qm);bm=dm*inv(qm);wm=fmm*inv(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)=xn(j)A(i-1);if j=0 | i=1cn(j,i)=0;elsecn(j,i)=(i-1)*xn(j)A(i-2);end

8、if j=0 | i=1 | i=2dn(j,i)=0;elsedn(j,i)=(i-2)*(i-1)*xn(j)A(i-3);endendfnn( j)=1/j;endan=cn*inv(qn);bn=dn*inv(qn);wn=fnn*inv(qn);为形状因子,m为配置点数,fm为权函数%正交多项式求根（适用于对称问题）function p=symm(a,m,fm) %afor 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);%权函数为 1elsec2=c2*(m+a/2+j+1);%权函数为

9、 1-xA2endc3=c3*(a/2+j);c4=c4*(1+j);endp(m+1-i)=c1*c2/c4/c3;endp(m+1)=1;% 为多项式系数向量，求出根后对对称问题还应开方才是零点p=sqrt(roots(p);%正交多项式求根(适用于非对称性问题)function p=unsymm(n,fn)if fn=0r(1)=(-1)An*n*(n+1);%权函数为 1elser(1)=(-1)An*n*(n+2);%权函数为 1-xendfor i=1:n-1if fn=0r(i+1)=(n-i)*(i+n+1)*r(i)/(i+1)/(i+1);%权函数为 1elser(i+1)

10、=(n-i)*(i+n+2)*r(i)/(i+1)/(i+1);%权函数为 1-xendendfor j=1:nP(n+1-j)=(-1)A(j+1)*r(j);endp(n+1)=(-1F(n+1);p=roots(p);应用正交配置法求解以下等温球形催化剂颗粒内反应物浓度分布，其浓度分布的数学模型为:1 dr2 drr r2 dC r drdC 0,dr1,C Cs36CR2解:(1)标准化令x r/R, y C/Cs代入微分方程及边界条件得：d 2 dy7 x 36yx2 dxdxx 0,曳 0dxx 1, y 1(2)离散化N 1Bji yi 36yj 0 j 1, 2, N 1i 1

11、(3)转化为代数方程组(以 N 3为例)B1136Bi2B13B14必0B21B22 36B23B24V20B31B32B3336B34y30B41B42B43B44 36V40因为yN 1 V41 ,所以整理上式得：B1136B12B13B14B21B22 36B23y1B24B31B32B3336y2B34B41B42B43y3B4436本例中的代数方程组为线性方程组，可采用线性方程组的求解方法；若为非线性方程组则采用相应的方法求解。命令：N=3,权函数为1-x 2am,bm,wm,an,bn,wn=collmatrix(3,3,1,3,1);(只用对称性配置矩阵)b1=bm;for i=

12、1:4b1(i,i)=bm(i,i)-36;enda0=b1(1:4,1:3);b0=-b1(1:4,4);y=a0b0;y(4)=1;p=exam31(3,3);(注意要对文件修改权函数为1-x2)x=0.3631,0.6772,0.8998,1;%零点plot(x,y,o)hold onx0=0:0.1:1;% 真实解y0=sinh(6*x0)./x0/sinh(6);plot(x0,y0,r)若权函数改为1,则以下语句修改，其他不变am,bm,wm,an,bn,wn=collmatrix(3,3,0,3,1);(只用对称性配置矩阵)p=exam31(3,3);(注意要对文件修改权函数为1

13、)x=0.4058,0.7415,0.9491,1;% 零点计算结果:权函数为1- x 2U正交彘量法Il -i权函数为11正交配置法 真实解0.90.80.70.6y 0.50.40.30.2yi V2V2 4yi2yi 0i, yi i e边值问题的MatLab 解法y 4y 0 x 12y 0 1, y 1 e精确解:2xy e函数：(collfunl.m )function f=collfun1(x,y)f(i)=y(2);f(2)=4*y(1);f=f(1)；f(2)；(collbcl.m )function f=collbc1(a,b)f=a(1)-1;b(1)-exp(2);命令

14、：solinit=bvpinit(0:0.1:1,1,1)sol=bvp4c(collfun1,collbc1,solinit)plot(sol.x,sol.y)hold onplot(sol.x,exp(2*sol.x),*)真实解15 I.yi Iy2真实解2 xy x 1 y 2y 1 x e 0 x 1y 02, y 11/eyi y22cxy21 x e 2y1 x 1 y2y1 02, y1 11/e精确解:函数：(collfun2.m ) function f=collfun2(x,y)f(i)=y(2);f(2)=(1-x.A2).*exp(-x)+2*y(1)-(x+1).*

15、y(2);f=f(1)；f(2)；(collbc2.m )function f=collbc2(a,b)f=a(2)-2;b(2)-exp(-1);命令：solinit=bvpinit(0:0.1:1,1,1);sol=bvp4c(collfun2,collbc2,solinit);plot(sol.x,sol.y)hold onplot(sol.x,(sol.x-1).*exp(-sol.x),*)真实解2 y 1y y - - 2 In x 1 x 2x x2y 1 y 10, y 23/2精确解：y22y1 1y x In xy1y21 2 In x xy210,y1xy2 23/2函数

16、：(collfun3.m )function f=collfun3(x,y)f(1)=y(2);f(2)=(2-log(x)./x+y(1)./x-y(2).A2;f=f;f(2);(collbc3.m )function f=collbc3(a,b)f=2*a(1)-a(2);b(2)-1.5;命令：solinit=bvpinit(1:0.1:2,1,1);sol=bvp4c(collfun3,collbc3,solinit);plot(sol.x,sol.y)hold onplot(sol.x,sol.x+log(sol.x),*)真实解在260 C的基础面上，为促进传热在此表面上增加纯铝

17、的圆柱形肋片，其直径为25mm ,高为150mm ;该柱表面受到 16 C气流的冷却，气流 与肋片表面的对流传热系数为15W/m2 K ,肋端绝热；肋片的导热系数为236 W/m K,假设肋片的导热热阻与肋片表面的对流传热热阻相比可以忽略；试求肋片中的温度分布，及单个肋 片的散热量为多少?解：根据以上条件可知：导热热阻与对流热阻相比可以忽略，则在肋片径向上没有温 度分布，在轴向上存在温度分布，外界气流与肋片的对流传热则可转化为内热源；故该问题为导热系数为常数的一维稳定热传导，其导热微分方程为:边界条件为:d2tdx2hp t tACx 0 时，to 260 C （肋根）；x H 时，一 dxx

18、 H0 （肋端绝热）ch m x Hhp 白 八分析斛：t tt0 t m J；传热重： Q0 ch mH. ACACx0这是两点边值的常微分方程求解问题，故可转化为如下形式:V2y1 0260,yiy2hp y yACy2 H0, 0 x H函数：(leipianfun.m leipianbc.m )%圆柱形肋片(常微分方程组)function f=leipianfun(x,y)f(i)=y(2);f(2)=15*pi*0.025/236/pi/0.025A2*4*(y(1)-16);f=f(1)；f(2)；%圆柱形肋片（边界条件）function f=leipianbc(a,b) f=a(

19、1)-260;b(2);命令：solinit=bvpinit(0:0.01:0.150,1,1);sol=bvp4c(leipianfun,leipianbc,solinit);%sol.y 中每行对应 sol.x 节点的因变量值即：第一行为y1 ,第二行为y2值，依此类推；故第一行为函数值，第二行对应的一阶导数。plot(sol.x,sol.y(1,:)%以下为分析解m=sqrt(15*pi*0.025/236/(pi/4*0.025A2)c1=(exp(m*(sol.x-0.15)+exp(-m*(sol.x-0.15)/2;c2=(exp(m*(0.15)+exp(-m*(0.15)/2

20、;t=16+(260-16)*c1/c2;hold onplot(sol.x,t,*)%计算传热量q=-236*pi*0.025A2/4*sol.y(2,1);计算结果:q=40.1052W在直彳至为20mm 的圆管外安装环形肋片，其表面温度为 260 C ,肋片导热系数为245 W/m K ,置于16 C、对流传热系数为150 W/m K的气流中；试根据单个环肋的传热量大小确定适宜的肋片高度和肋片厚度；并给出肋高为0.01m ,肋厚为0.0003m 环肋 的温度分布。解：近似为：导热系数为常数的一维稳定热传导，其导热微分方程为:1 d dt hp t t 2h t tr r dr drAC边

21、界条件为：r %时，t0 260 C （肋根）；r 2 r H时，一 0 （肋端绝热）这是两点边值的常微分方程求解问题，故可转化为如下形式：必 V22h y y y2y2xy1 rI260, V2 r20, r x 万函数：(huanleifun.mhuanleibc.m )%环形肋片(常微分方程组)function f=huanleifun(x,y,h,nada,delta,t0,tf)f(i)=y(2);f(2)=2*h/nada/delta*(y(1)-tf)-y(2)./x;f=f(1)；f(2)；%环形肋片(边界条件)function f=huanleibc(a,b,h,nada,d

22、elta,t0,tf) f=a(1)-t0;b(2);命令：(hq.m )%环肋不同肋高对散热量的影响function q,sol=hqh=150;nada=45;delta=0.0003;t0=260;tf=16;r=0.01;H=0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01,0.012,0.014,0.016,0.018,0.02,0.03;x=r+H;for i=1:length(x)solinit=bvpinit(r:0.001:x(i),1,1);sol=bvp4c(huanleifun,huanleibc,so

23、linit,口,h,nada,delta,t0,tf);q(i)=-nada*2*pi*r*delta*sol.y(2,1);endH=0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01,0.012,0.014,0.016,0.018,0.02,0.03;q,sol=hq;plot(H,q)计算结果：(肋厚为0.3mm )由图可知，适宜的肋高可取0.0050.015m。命令：（dq.m ）%环肋不同肋厚对散热量的影响function q,sol=dq h=150;nada=45;delta=0.0001,0.0002,0.000

24、3,0.0004,0.0005,0.0006,0.0008,0.0009,0.001,0.0015,0.002,0.0025,0.003,0.0035,0.004,0.005;t0=260;tf=16;r=0.01;x=r+0.01;for i=1:length(delta)solinit=bvpinit(r:0.001:x,1,1);sol=bvp4c(huanleifun,huanleibc,solinit,口,h,nada,delta(i),t0,tf);q(i)=-nada*2*pi*r*delta(i)*sol.y(2,1);enddelta=0.0001,0.0002,0.0003,0.0004,0.0005,0.0006,0.0008,0.0009,0.001,0.0015,0.002,0.0025,0.003,0.0035,0.004,0.005;q,sol=dq;plot(delta,q)计算结果：（肋高为10mm ）由图可知，适宜的肋厚可取0.00050.002m ,而在同一基础面上肋厚越大，则肋片数目越少，虽然肋厚增