常微分方程组边值.docx

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; 应用打靶法求解下列边值问题：

2、 解：将其转化为常微分方程组的初值问题命令：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)函数：（odebvp.m）%边值常微分方程(组)函数function

3、 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） 同上例：若划分为10个区间，则：函数：（difference.m）%有限差分法求常微分方程的边值问

4、题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);计算结果：x0=0:0.1:10;y0=32*(cos(5

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

6、rix(a,m,fm,n,fn)x0=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);

7、wm=fmm*inv(qm);x1=unsymm(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)

8、=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: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

9、+1)=1;%为多项式系数向量,求出根后对对称问题还应开方才是零点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)

10、=(-1)(j+1)*r(j);endp(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;e

11、nda0=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）x=0.4058,0.7415,0.949

12、1,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),*) 真实解 精确解：函数

13、：（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),*) 真实解 精确解：函数：（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