常微分方程组(边值)

,.常微分方程组边值问题解法打靶法ShootingMethod（shooting.m）谢谢阅读%打靶法求常微分方程的边值问题function[x,a,b,n]=shooting(fun,x0,xn,eps)谢谢阅读ifnargin<3eps=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&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;应用打靶法求解下列边值问题：d2yy8dx24y00y100解：将其转化为常微分方程组的初值问题dyy1dx2dyy2814dxy0y10dytdxx0命令：x0=[0:0.1:10];,.y0=32*((cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); 真实解感谢阅读plot(x0,y0,'r')holdon[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)),.函数：（odebvp.m）%边值常微分方程(组)函数functionf=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)感谢阅读,.计算结果：（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);精品文档放心下载holdonplot(x0,y0,’o’),.有限差分法FiniteDifferenceMethods FDM（difference.m）谢谢阅读同上例：d2yyy2yyi18y8i1iidx24h24y2h24yy8h2i1i1i若划分为10个区间，则：,.2h214h2y121204y18h28h2h2y8h2121yn28h204n112h24函数：（difference.m）%有限差分法求常微分方程的边值问题function[x,y]=difference(x0,xn,y0,yn,n)精品文档放心下载h=(xn-x0)/n;a=eye(n-1)*(-(2-h^2/4));fori=1:n-2a(i,i+1)=1;a(i+1,i)=1;endb=ones(n-1,1)*8*h^2;b(1)=b(1)-0;b(n-1)=b(n-1)-0;yy=a\b;,.x(1)=x0;y(1)=y0;fori=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')holdon[x,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,’-.’),.正交配置法OrthogonalCollocatioinMethods CM精品文档放心下载构造正交矩阵函数（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-x^2)fori=1:mxm(i)=x0(m+1-i);endxm(m+1)=1;forj=1:m+1fori=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));感谢阅读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;fori=2:n+1xn(i)=x1(n+2-i);endxn(n+2)=1;forj=1:n+2fori=1:n+2qn(j,i)=xn(j)^(i-1);ifj==0|i==1,.cn(j,i)=0;elsecn(j,i)=(i-1)*xn(j)^(i-2);谢谢阅读endifj==0|i==1|i==2dn(j,i)=0;elsedn(j,i)=(i-2)*(i-1)*xn(j)^(i-3);精品文档放心下载endendfnn(j)=1/j;endan=cn*inv(qn);bn=dn*inv(qn);wn=fnn*inv(qn);%正交多项式求根(适用于对称问题)functionp=symm(a,m,fm)%a为形状因子,m为配置点数,fm为权函数精品文档放心下载,.fori=1:mc1=1;c2=1;c3=1;c4=1;forj=0:i-1c1=c1*(-m+j);iffm==0c2=c2*(m+a/2+j);%权函数为1elsec2=c2*(m+a/2+j+1);%权函数为1-x^2感谢阅读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));,.%正交多项式求根(适用于非对称性问题)functionp=unsymm(n,fn)iffn==0r(1)=(-1)^n*n*(n+1);%权函数为1感谢阅读elser(1)=(-1)^n*n*(n+2);%权函数为1-x谢谢阅读endfori=1:n-1iffn==0r(i+1)=(n-i)*(i+n+1)*r(i)/(i+1)/(i+1);%权函数为1精品文档放心下载elser(i+1)=(n-i)*(i+n+2)*r(i)/(i+1)/(i+1);%权函数为1-x感谢阅读endendforj=1:np(n+1-j)=(-1)^(j+1)*r(j);end,.p(n+1)=(-1)^(n+1);p=roots(p);应用正交配置法求解以下等温球形催化剂颗粒内反应物浓度分布，其浓度分布的数学谢谢阅读模型为：1ddC36Cr2r2drdrR2dCr0,0drr1,CCS解：（1）标准化令xr/R，yC/CS代入微分方程及边界条件得：谢谢阅读1ddy36yx2x2dxdxdyx0,0dxx1,y1（2）离散化N1j1,2,N1By36y0jiij1（3）转化为代数方程组（以N3为例）,.B36BBBy011123613141B21BBB24y20B22B23By0B36B31B3233B34304142B36y44344因为yy1，所以整理上式得：N14B36BByB11B12361314BB1BB2122B2336yB2431B32234BB33yB3B3641424344本例中的代数方程组为线性方程组，可采用线性方程组的求解方法；若为非线性方程组则谢谢阅读采用相应的方法求解。命令：N=3,权函数为1-x2[am,bm,wm,an,bn,wn]=collmatrix(3,3,1,3,1);（只用对称性配置矩阵）感谢阅读b1=bm;fori=1:4b1(i,i)=bm(i,i)-36;enda0=b1(1:4,1:3);b0=-b1(1:4,4);y=a0\b0;y(4)=1;p=exam31(3,3);（注意要对文件修改权函数为1-x2）感谢阅读,.x=[0.3631,0.6772,0.8998,1]; %零点精品文档放心下载plot(x,y,'o')holdonx0=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.9491,1]; %零点谢谢阅读计算结果：权函数为1-x2,.权函数为11正交配置法0.9真实解y 0.2,.边值问题的MatLab解法y4y0x1yy12e2y4y21y01,y11,ye20111,.精确解：y e2x函数：（collfun1.m）functionf=collfun1(x,y)f(1)=y(2);f(2)=4*y(1);f=[f(1);f(2)];（collbc1.m）functionf=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)holdonplot(sol.x,exp(2*sol.x),'*') 真实解精品文档放心下载,.2x0x1y02,y11/e精确解：yx1ex

y 2

y1y2x1y1x2ex2y1211,.函数：（collfun2.m）functionf=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）functionf=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)holdonplot(sol.x,(sol.x-1).*exp(-sol.x),'*') 真实解感谢阅读,.yy22y1精确解：

y1yy1122lnx1x2yyy2xx2lnx12xx2y10,y23/22y1y10,y23/2122yxlnx函数：（collfun3.m）functionf=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）functionf=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)holdonplot(sol.x,sol.x+log(sol.x),'*') 真实解感谢阅读在260C的基础面上，为促进传热在此表面上

,.H260C增加纯铝的圆柱形肋片，其直径为25mm，高为150mm；该柱表面受到16C气流的冷却，气流谢谢阅读与肋片表面的对流传热系数为15W/m2K，肋端谢谢阅读16C绝热；肋片的导热系数为236W/mK，假设肋片 气流感谢阅读的导热热阻与肋片表面的对流传热热阻相比可以忽略；试求肋片中的温度分布，及单个肋精品文档放心下载片的散热量为多少？,.解：根据以上条件可知：导热热阻与对流热阻相比可以忽略，则在肋片径向上没有温精品文档放心下载度分布，在轴向上存在温度分布，外界气流与肋片的对流传热则可转化为内热源；故该问谢谢阅读题为导热系数为常数的一维稳定热传导，其导热微分方程为：谢谢阅读2tdhpttdx2AC边界条件为：x0时，t0260C（肋根）；xH时，dt0（肋端绝热）。dxxH分析解：ttttchmxH，mhp；传热量：QAdt0chmHACCdxx0这是两点边值的常微分方程求解问题，故可转化为如下形式：yy12hpyyy2A260,yC0xHy0H0,12函数：（leipianfun.m leipianbc.m）谢谢阅读%圆柱形肋片(常微分方程组)functionf=leipianfun(x,y)感谢阅读f(1)=y(2);f(2)=15*pi*0.025/236/pi/0.025^2*4*(y(1)-16);精品文档放心下载f=[f(1);f(2)];%圆柱形肋片(边界条件),.functionf=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节点的因变感谢阅读量值即：第一行为y1，第二行为y2值，依此类推；故第一行为函数值，第二行对应的一精品文档放心下载阶导数。plot(sol.x,sol.y(1,:))%以下为分析解m=sqrt(15*pi*0.025/236/(pi/4*0.025^2))精品文档放心下载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;holdonplot(sol.x,t,'*')%计算传热量q=-236*pi*0.025^2/4*sol.y(2,1);谢谢阅读,.计算结果：q=40.1052W,.在直径为20mm的圆管外安装环形肋片，其表面温度为260C，肋片导热系数为感谢阅读45W/mK，置于16C、对流传热系数为150W/m2K的气流中；试根据单个环肋的传谢谢阅读热量大小确定适宜的肋片高度和肋片厚度；并给出肋高为0.01m，肋厚为0.0003m环肋精品文档放心下载的温度分布。解：近似为：导热系数为常数的一维稳定热传导，其导热微分方程为：谢谢阅读1ddthptt2httrArdrdrC边界条件为：rr时，t260C（肋根）；rrrH时，dt0（肋端绝热）。1021drrr2这是两点边值的常微分方程求解问题，故可转化为如下形式：yy12yy2hyy22xyr260,yr0,rxr112212函数：（huanleifun.m huanleibc.m）精品文档放心下载,.%环形肋片(常微分方程组)functionf=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)];%环形肋片(边界条件)functionf=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;fori=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.005~0.015m。谢谢阅读,.命令：（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.0精品文档放心下载02,0.0025,0.003,0.0035,0.004,0.005];感谢阅读t0=260;tf=16;r=0.01;x=r+0.01;fori=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.0谢谢阅读02,0.0025,0.003,0.0035,0.004,0.005];[q,sol]=dq;plot(delta,q),.计算结果：（肋高为10mm）由图可知，适宜的肋厚可取0.0005~0.002m，而在同一基础面上肋厚越大，则肋片感谢阅读数目越少，虽然肋厚增加散