常微分方程组边值问题解法打靶法Shooting Method （shooting.m ）%打靶法求常微分方程的边值问题function [x,a,b,n]=shooting(fun,x0,xn,eps)if nargin<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;应用打靶法求解下列边值问题：()()⎪⎪⎩⎪⎪⎨⎧==-=010004822y y y dxy d解：将其转化为常微分方程组的初值问题()⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧==-==t dx dy y y y dx dy y dx dy x 0011221048命令：x0=[0:0.1:10];y0=32*((cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); 真实解plot(x0,y0,'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))函数：（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]=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 ）同上例：4822y dx y d -=⇒482211i i i i y h y y y -=+--+ 2121842h y y h y i i i =+⎪⎪⎭⎫ ⎝⎛--+- 若划分为10个区间，则： ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛--⎪⎪⎭⎫ ⎝⎛--⎪⎪⎭⎫ ⎝⎛--⎪⎪⎭⎫ ⎝⎛----08880842114211421142222212212222h h h h y y y y h h h h n n M M O O O函数：（difference.m ）%有限差分法求常微分方程的边值问题function [x,y]=difference(x0,xn,y0,yn,n)h=(xn-x0)/n;a=eye(n-1)*(-(2-h^2/4));for i=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;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 on[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,’-.’)正交配置法Orthogonal Collocatioin Methods 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)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)^(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;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)^(i-1);if j==0 | i==1cn(j,i)=0;elsecn(j,i)=(i-1)*xn(j)^(i-2);endif j==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);%正交多项式求根(适用于对称问题)function p=symm(a,m,fm) %a为形状因子,m为配置点数,fm为权函数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);%权函数为1elsec2=c2*(m+a/2+j+1);%权函数为1-x^2endc3=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)^n*n*(n+1);%权函数为1elser(1)=(-1)^n*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)=(n-i)*(i+n+2)*r(i)/(i+1)/(i+1);%权函数为1-xendendfor j=1:np(n+1-j)=(-1)^(j+1)*r(j);endp(n+1)=(-1)^(n+1);p=roots(p);应用正交配置法求解以下等温球形催化剂颗粒内反应物浓度分布，其浓度分布的数学模型为：⎪⎪⎪⎩⎪⎪⎪⎨⎧=====⎪⎭⎫ ⎝⎛S C C r dr dC r R C dr dC r dr d r ,10,0361222 解：（1）标准化令R r x /=，S C C y /=代入微分方程及边界条件得：⎪⎪⎪⎩⎪⎪⎪⎨⎧=====⎪⎭⎫ ⎝⎛1,10,036122y x dx dy x y dx dy x dx d x（2）离散化03611=-∑+=j N i i ji y y B1,2,1+=N j Λ（3）转化为代数方程组（以3=N 为例）⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----000036363636432144434241343332312423222114131211y y y y B B B B B B B B B B B B B B B B 因为141==+y y N ，所以整理上式得：⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡+----=⎥⎥⎦⎤⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---3636363644342414321434241333231232221131211B B B B y y y B B B B B B B B B B B B 本例中的代数方程组为线性方程组，可采用线性方程组的求解方法；若为非线性方程组则采用相应的方法求解。

命令：N=3,权函数为1-x2 [am,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=-b1(1:4,4);y=a0\b0;y(4)=1;p=exam31(3,3);（注意要对文件修改权函数为1-x 2）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.9491,1]; %零点计算结果：权函数为1- x2权函数为1y 正交配置法真实解边值问题的MatLab 解法()()⎩⎨⎧==<<=''21,10104e y y x y y ⇒ ()()⎪⎩⎪⎨⎧==='='21112211,104ey y y y y y 精确解：x e y 2=函数：（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),'*') 真实解()()()()⎩⎨⎧='='<<-=-'++''-e y y x e x y y x y x /11,20101212 ⇒ ()()()()⎪⎩⎪⎨⎧='='+-+-='='-e y y y x y e x y y y x /11,2012111212221 精确解：()x e x y --=1函数：（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),'*') 真实解()()()()⎪⎩⎪⎨⎧='='-<<-=-'+''2/32,011221ln 212y y y x x xx y y y ⇒ ()()()()⎪⎪⎩⎪⎪⎨⎧==--+-='='2/32,0112ln 21221221221y y y y x y x x y y y 精确解：x x y ln +=函数：（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),'*') 真实解在260C ︒的基础面上，为促进传热在此表面上增加纯铝的圆柱形肋片，其直径为25mm ，高为150mm ；该柱表面受到16C ︒气流的冷却，气流与肋片表面的对流传热系数为15K m W ⋅2/，肋端绝热；肋片的导热系数为236K m W ⋅/，假设肋片的导热热阻与肋片表面的对流传热热阻相比可以忽略；试求肋片中的温度分布，及单个肋片的散热量为多少？解：根据以上条件可知：导热热阻与对流热阻相比可以忽略，则在肋片径向上没有温度分布，在轴向上存在温度分布，外界气流与肋片的对流传热则可转化为内热源；故该问题为导热系数为常数的一维稳定热传导，其导热微分方程为：()C A t t hp dxt d λλ∞-=Φ=&22 边界条件为：0=x 时，C 2600︒=t （肋根）；H x =时，0==Hx dxdt（肋端绝热）。