常微分方程組(邊值).doc

常微分方程組邊值問題解法打靶法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; 應(yīng)用打靶法求解下列邊值問題： 解：將其轉(zhuǎn)化為常微分方程組的初值問題命令：x0=0:0.1:10;y0=32*(cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); 真實(shí)解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)函數(shù)：（odebvp.m）%邊值常微分方程(組)函數(shù)function f=odebvp(x,y)f(1)=y(2);f(2)=8-y(1)/4;f=f(1);f(2);命令：t,x,y,n=shooting(odebvp,10,0,1e-3)計(jì)算結(jié)果：（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個(gè)區(qū)間，則：函數(shù)：（difference.m）%有限差分法求常微分方程的邊值問題function x,y=difference(x0,xn,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);計(jì)算結(jié)果：x0=0:0.1:10;y0=32*(cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); 真實(shí)解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構(gòu)造正交矩陣函數(shù)（collmatrix.m）%正交配置矩陣(均用矩陣法求對(duì)稱性與非對(duì)稱性正交配置矩陣)function am,bm,wm,an,bn,wn=collmatrix(a,m,fm,n,fn)x0=symm(a,m,fm); %a為形狀因子;m為零點(diǎn)數(shù);fm為對(duì)稱的權(quán)函數(shù)(0為權(quán)函數(shù)1,非0為權(quán)函數(shù)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=unsymm(n,fn); %n為零點(diǎn)數(shù);fn為非對(duì)稱的權(quán)函數(shù)(0為權(quán)函數(shù)1,非0為權(quán)函數(shù)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(qn);bn=dn*inv(qn);wn=fnn*inv(qn);%正交多項(xiàng)式求根(適用于對(duì)稱問題)function p=symm(a,m,fm) %a為形狀因子,m為配置點(diǎn)數(shù),fm為權(quán)函數(shù)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);%權(quán)函數(shù)為1 else c2=c2*(m+a/2+j+1);%權(quán)函數(shù)為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;%為多項(xiàng)式系數(shù)向量,求出根后對(duì)對(duì)稱問題還應(yīng)開方才是零點(diǎn)p=sqrt(roots(p);%正交多項(xiàng)式求根(適用于非對(duì)稱性問題)function p=unsymm(n,fn)if fn=0 r(1)=(-1)n*n*(n+1);%權(quán)函數(shù)為1else r(1)=(-1)n*n*(n+2);%權(quán)函數(shù)為1-xendfor i=1:n-1 if fn=0 r(i+1)=(n-i)*(i+n+1)*r(i)/(i+1)/(i+1);%權(quán)函數(shù)為1 else r(i+1)=(n-i)*(i+n+2)*r(i)/(i+1)/(i+1);%權(quán)函數(shù)為1-x endendfor j=1:n p(n+1-j)=(-1)(j+1)*r(j);endp(n+1)=(-1)(n+1);p=roots(p); 應(yīng)用正交配置法求解以下等溫球形催化劑顆粒內(nèi)反應(yīng)物濃度分布，其濃度分布的數(shù)學(xué)模型為： 解： （1）標(biāo)準(zhǔn)化 令，代入微分方程及邊界條件得： （2）離散化 （3）轉(zhuǎn)化為代數(shù)方程組（以為例）因?yàn)?，所以整理上式得：本例中的代?shù)方程組為線性方程組，可采用線性方程組的求解方法；若為非線性方程組則采用相應(yīng)的方法求解。命令：N=3,權(quán)函數(shù)為1-x2am,bm,wm,an,bn,wn=collmatrix(3,3,1,3,1);（只用對(duì)稱性配置矩陣）b1=bm;for i=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);（注意要對(duì)文件修改權(quán)函數(shù)為1-x2）x=0.3631,0.6772,0.8998,1; %零點(diǎn)plot(x,y,o)hold onx0=0:0.1:1; %真實(shí)解y0=sinh(6*x0)./x0/sinh(6);plot(x0,y0,r)若權(quán)函數(shù)改為1，則以下語句修改，其他不變am,bm,wm,an,bn,wn=collmatrix(3,3,0,3,1);（只用對(duì)稱性配置矩陣）p=exam31(3,3);（注意要對(duì)文件修改權(quán)函數(shù)為1）x=0.4058,0.7415,0.9491,1; %零點(diǎn)計(jì)算結(jié)果：權(quán)函數(shù)為1- x2權(quán)函數(shù)為1邊值問題的MatLab解法 精確解：函數(shù)：（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),*) 真實(shí)解 精確解：函數(shù)：（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),*) 真實(shí)解 精確解：函數(shù)：（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),*) 真實(shí)解氣流 在260的基礎(chǔ)面上，為促進(jìn)傳熱在此表面上增加純鋁的圓柱形肋片，其直徑為25mm，高為150mm；該柱表面受到16氣流的冷卻，氣流與肋片表面的對(duì)流傳熱系數(shù)為15，肋端絕熱；肋片的導(dǎo)熱系數(shù)為236，假設(shè)肋片的導(dǎo)熱熱阻與肋片表面的對(duì)流傳熱熱阻相比可以忽略；試求肋片中的溫度分布，及單個(gè)肋片的散熱量為多少？ 解：根據(jù)以上條件可知：導(dǎo)熱熱阻與對(duì)流熱阻相比可以忽略，則在肋片徑向上沒有溫度分布，在軸向上存在溫度分布，外界氣流與肋片的對(duì)流傳熱則可轉(zhuǎn)化為內(nèi)熱源；故該問題為導(dǎo)熱系數(shù)為常數(shù)的一維穩(wěn)定熱傳導(dǎo)，其導(dǎo)熱微分方程為：邊界條件為：時(shí)，（肋根）；時(shí)，（肋端絕熱）。 分析解：，；傳熱量： 這是兩點(diǎn)邊值的常微分方程求解問題，故可轉(zhuǎn)化為如下形式：函數(shù)：（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中每行對(duì)應(yīng)sol.x節(jié)點(diǎn)的因變量值即：第一行為y1，第二行為y2值，依此類推；故第一行為函數(shù)值，第二行對(duì)應(yīng)的一階導(dǎo)數(shù)。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,*)%計(jì)算傳熱量q=-236*pi*0.0252/4*sol.y(2,1);計(jì)算結(jié)果：q=40.1052W 在直徑為20mm的圓管外安裝環(huán)形肋片，其表面溫度為260，肋片導(dǎo)熱系數(shù)為45，置于16、對(duì)流傳熱系數(shù)為150的氣流中；試根據(jù)單個(gè)環(huán)肋的傳熱量大小確定適宜的肋片高度和肋片厚度；并給出肋高為0.01m，肋厚為0.0003m環(huán)肋的溫度分布。 解：近似為：導(dǎo)熱系數(shù)為常數(shù)的一維穩(wěn)定熱傳導(dǎo)，其導(dǎo)熱微分方程為：邊界條件為：時(shí)，（肋根）；時(shí)，（肋端絕熱）。 這是兩點(diǎn)邊值的常微分方程求解問題，故可轉(zhuǎn)化為如下形式：函數(shù)：（huanleifun.m huanleibc.m）%環(huán)形肋片(常微分方程組)function f=huanleifun(x,y,h,nada,delta,t0,tf)f(1)=y(2);f(2)=2*h/nada/delta*(y(1)-tf)-y(2)./x;f=f(1);f(2);%環(huán)形肋片(邊界條件)function f=huanleibc(a,b,h,nada,delta,t0,tf)f=a(1)-t0;b(2);命令：（hq.m）%環(huán)肋不同肋高對(duì)散熱量的影響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,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)計(jì)算結(jié)果：（肋厚為0.3mm） 由圖可知，適宜的肋高可取0.0050.015m。命令：（dq.m）%環(huán)肋不同肋厚對(duì)散熱量的影響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)計(jì)算結(jié)果：（肋高為10mm） 由圖可知，適宜的肋厚可取0.00050.002m，而在同一基礎(chǔ)面上肋厚越大，則肋片數(shù)目越