常微分方程组边值

1、常微分方程组边值问题解法打靶法Shooting Method（shooting.m）%打靶法求常微分方程的边值问题function x,a,b,n=shooting(fun,x0,xn,eps)if nargin<3 eps=1e-3;endx1=x0+rand;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=1;while (norm(c1-xn)>=eps & no

2、rm(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; 应用打靶法求解下列边值问题： 解：将其转化为常微分方程组的初值问题命令：x0=0:0.1:10;y0=32*(cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); 真实解plot(x0,y0,'r')h

3、old 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)函数：（odebvp.m）%边值常微分方程(组)函数function f=odebvp(x,y)f(1)=y(2);f(2)=8-y(1)/4;f

4、=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)有限差分法Finite Difference Methods FDM（difference.m） 同上例：若划分为10个区间，则：函数：（difference.m）%有限差分法求常微分方程的边值问题function x,y=difference(x0,x

5、n,y0,yn,n)h=(xn-x0)/n;a=eye(n-1)*(-(2-h2/4);for i=1:n-2 a(i,i+1)=1; a(i+1,i)=1;endb=ones(n-1,1)*8*h2;b(1)=b(1)-0;b(n-1)=b(n-1)-0;yy=ab;x(1)=x0;y(1)=y0;for i=2:n x(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/

6、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);plot(x,y,-.)正交配置法Orthogonal Collocatioin Methods CM构造正交矩阵函数（collmatrix.m）%正交配置矩阵(均用矩阵法求对称性与非对称性正交配置矩阵)function am,bm,wm,an,bn,wn=collmatrix(a,m,fm,n,fn)x0=

7、symm(a,m,fm); %a为形状因子;m为零点数;fm为对称的权函数(0为权函数1,非0为权函数1-x2)for i=1:m xm(i)=x0(m+1-i);endxm(m+1)=1;for j=1:m+1 for i=1:m+1 qm(j,i)=xm(j)(2*i-2); cm(j,i)=(2*i-2)*xm(j)(2*i-3); dm(j,i)=(2*i-2)*(2*i-3+(a-1)*xm(j)(2*i-3+(a-1)-1-(a-1); end fmm(j)=1/(2*j-2+a);endam=cm*inv(qm);bm=dm*inv(qm);wm=fmm*inv(qm);x1=u

8、nsymm(n,fn); %n为零点数;fn为非对称的权函数(0为权函数1,非0为权函数1-x)xn(1)=0;for i=2:n+1 xn(i)=x1(n+2-i);endxn(n+2)=1;for j=1:n+2 for i=1:n+2 qn(j,i)=xn(j)(i-1); if j=0 | i=1 cn(j,i)=0; else cn(j,i)=(i-1)*xn(j)(i-2); end if j=0 | i=1 | i=2 dn(j,i)=0; else dn(j,i)=(i-2)*(i-1)*xn(j)(i-3); end end fnn(j)=1/j;endan=cn*inv(q

9、n);bn=dn*inv(qn);wn=fnn*inv(qn);%正交多项式求根(适用于对称问题)function p=symm(a,m,fm) %a为形状因子,m为配置点数,fm为权函数for i=1:m c1=1; c2=1; c3=1; c4=1; for j=0:i-1 c1=c1*(-m+j); if fm=0 c2=c2*(m+a/2+j);%权函数为1 else c2=c2*(m+a/2+j+1);%权函数为1-x2 end c3=c3*(a/2+j); c4=c4*(1+j); end p(m+1-i)=c1*c2/c4/c3;endp(m+1)=1;%为多项式系数向量,求出根

10、后对对称问题还应开方才是零点p=sqrt(roots(p);%正交多项式求根(适用于非对称性问题)function p=unsymm(n,fn)if fn=0 r(1)=(-1)n*n*(n+1);%权函数为1else r(1)=(-1)n*n*(n+2);%权函数为1-xendfor i=1:n-1 if fn=0 r(i+1)=(n-i)*(i+n+1)*r(i)/(i+1)/(i+1);%权函数为1 else r(i+1)=(n-i)*(i+n+2)*r(i)/(i+1)/(i+1);%权函数为1-x endendfor j=1:n p(n+1-j)=(-1)(j+1)*r(j);end

11、p(n+1)=(-1)(n+1);p=roots(p); 应用正交配置法求解以下等温球形催化剂颗粒内反应物浓度分布，其浓度分布的数学模型为： 解： （1）标准化 令，代入微分方程及边界条件得： （2）离散化 （3）转化为代数方程组（以为例）因为，所以整理上式得：本例中的代数方程组为线性方程组，可采用线性方程组的求解方法；若为非线性方程组则采用相应的方法求解。命令：N=3,权函数为1-x2am,bm,wm,an,bn,wn=collmatrix(3,3,1,3,1);（只用对称性配置矩阵）b1=bm;for i=1:4b1(i,i)=bm(i,i)-36;enda0=b1(1:4,1:3);b0

12、=-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）x=0.4058,0.7415,0.94

13、91,1; %零点计算结果：权函数为1- x2权函数为1边值问题的MatLab解法 精确解：函数：（collfun1.m）function f=collfun1(x,y)f(1)=y(2);f(2)=4*y(1);f=f(1);f(2);（collbc1.m）function f=collbc1(a,b)f=a(1)-1;b(1)-exp(2);命令：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),'*')

14、 真实解 精确解：函数：（collfun2.m）function f=collfun2(x,y)f(1)=y(2);f(2)=(1-x.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);命令：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),'*&

15、#39;) 真实解 精确解：函数：（collfun3.m）function f=collfun3(x,y)f(1)=y(2);f(2)=(2-log(x)./x+y(1)./x-y(2).2;f=f(1);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),'*') 真实

16、解气流 在260的基础面上，为促进传热在此表面上增加纯铝的圆柱形肋片，其直径为25mm，高为150mm；该柱表面受到16气流的冷却，气流与肋片表面的对流传热系数为15，肋端绝热；肋片的导热系数为236，假设肋片的导热热阻与肋片表面的对流传热热阻相比可以忽略；试求肋片中的温度分布，及单个肋片的散热量为多少？ 解：根据以上条件可知：导热热阻与对流热阻相比可以忽略，则在肋片径向上没有温度分布，在轴向上存在温度分布，外界气流与肋片的对流传热则可转化为内热源；故该问题为导热系数为常数的一维稳定热传导，其导热微分方程为：边界条件为：时，（肋根）；时，（肋端绝热）。 分析解：，；传热量： 这是两点边值的常微

17、分方程求解问题，故可转化为如下形式：函数：（leipianfun.m leipianbc.m）%圆柱形肋片(常微分方程组)function f=leipianfun(x,y)f(1)=y(2);f(2)=15*pi*0.025/236/pi/0.0252*4*(y(1)-16);f=f(1);f(2);%圆柱形肋片(边界条件)function f=leipianbc(a,b)f=a(1)-260;b(2);命令：solinit=bvpinit(0:0.01:0.150,1,1);sol=bvp4c(leipianfun,leipianbc,solinit); %sol.y中每行对应sol.x节

18、点的因变量值即：第一行为y1，第二行为y2值，依此类推；故第一行为函数值，第二行对应的一阶导数。plot(sol.x,sol.y(1,:)%以下为分析解m=sqrt(15*pi*0.025/236/(pi/4*0.0252)c1=(exp(m*(sol.x-0.15)+exp(-m*(sol.x-0.15)/2;c2=(exp(m*(0.15)+exp(-m*(0.15)/2;t=16+(260-16)*c1/c2;hold onplot(sol.x,t,'*')%计算传热量q=-236*pi*0.0252/4*sol.y(2,1);计算结果：q=40.1052W 在直径为20

19、mm的圆管外安装环形肋片，其表面温度为260，肋片导热系数为45，置于16、对流传热系数为150的气流中；试根据单个环肋的传热量大小确定适宜的肋片高度和肋片厚度；并给出肋高为0.01m，肋厚为0.0003m环肋的温度分布。 解：近似为：导热系数为常数的一维稳定热传导，其导热微分方程为：边界条件为：时，（肋根）；时，（肋端绝热）。 这是两点边值的常微分方程求解问题，故可转化为如下形式：函数：（huanleifun.m huanleibc.m）%环形肋片(常微分方程组)function f=huanleifun(x,y,h,nada,delta,t0,tf)f(1)=y(2);f(2)=2*h/n

20、ada/delta*(y(1)-tf)-y(2)./x;f=f(1);f(2);%环形肋片(边界条件)function f=huanleibc(a,b,h,nada,delta,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;fo

21、r i=1:length(x)solinit=bvpinit(r:0.001:x(i),1,1);sol=bvp4c(huanleifun,huanleibc,solinit,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。命

22、令：（dq.m）%环肋不同肋厚对散热量的影响function q,sol=dqh=150;nada=45;delta=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;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，而在