常微分方程组边值Word格式.docx

1、应用打靶法求解下列边值问题:y 10 0解：将其转化为常微分方程组的初值问题命令:xO=O:O.1:1O;真实解y0=32*（cos（5）-1）/si n（ 5）*si n（x0/2）-cos（x0/2）+1）; plot（xO,yO,r）hold on x,y=ode45（odebvp,0,10,0,2 plot（x,y（:,1）x,y=ode45（,0,10,0,5plot（x,y（:,0,10,0,8,0,10,0,10函数：（odebvp.m）泌值常微分方程（组）函数function f=odebvp（x,y）f（1）=y（2）；f（2）=8-y（1）/4;f=f（1）；f（2）；命令

2、：t,x,y, n=shoot in g（,10,0,1e-3）计算结果：（eps=0.001 ）t=11.9524x0=0:1:10;hold onplot（xO,yO, o（differe nce.m ）有限差分法 Finite Differenee Methods FDM同上例:Yi i 2Yi yi iYi i 2若划分为10个区间，则：h- Yi Yi 1 8h24h22 12Y1Y28h2Yn 2Yn 1.2h%有限差分法求常微分方程的边值问题 function x,y=differe nce（xO,x n,yO,yn,n） h=（x n-xO）/n;a=eye（ n-1）*（-（

3、2-hA2/4）;for i=1: n-2a（i,i+1）=1;a（i+1,i）=1;b=o nes（ n-1,1）*8*hA2;b（1）=b（1）-0;b（n-1）=b（ n-1）-0;yy=ab;x（1）=x0;y（1）=y0;for i=2: nx（i）=x0+（i-1）*h;y（i）=yy（i-i）;x（n）=xn;y（n）=yn;x,y=differe nce（0,10,0,0,100）;plot（xO,yO,x,y=differe nce（0,10,0,0,5）;plot（x,y,.x,y=differe nce（0,10,0,0,10）;-x,y=differe nce（0,10

4、,0,0,50）;-.正交配置法 Orthogonal Collocatioin Methods CM构造正交矩阵函数（collmatrix.m ）%E交配置矩阵（均用矩阵法求对称性与非对称性正交配置矩阵 ）function am,bm,wm,an,bn,wn=collmatrix（a,m,fm,n,fn）xO=symm（a,m,fm）; %a为形状因子;m为零点数;fm为对称的权函数（0为权函数1,非0为权函数1-xA2）mxm（i）=x0（m+1-i）;xm（m+1）=1;for j=1:m+1qm（j,i）=xm（j）A（2*i-2）;cm（j,i）=（2*i-2）*xm（j）A（2*i

5、-3）;dm（j,i）=（2*i-2）*（2*i-3+（a-1）*xm（j）A（2*i-3+（a-1）-1-（a-1）;fmm（j）=1/（2*j-2+a）;am=cm*i nv（qm）;bm=dm*i nv（qm）;wm=fmm*i nv（qm）;x1=unsymm（n,fn）; %n 为零点数;fn为非对称的权函数（0为权函数1,非0为权函数1-x） xn （1）=0; n+1xn（ i）=x1（ n+2-i）;xn（n+2）=1; n+2qn （j,i）=x n（j）A（i-1）;if j=0 | i=1cn （j,i）=0;elsecn （j,i）=（i-1）*x n（j）A（i-2）

6、;if j=0 | i=1 | i=2dn （j,i）=O;dn （j,i）=（i-2）*（i-1）*x n（jF（i-3）;fnn （j）=1/j;an=cn *i nv（qn）;bn=dn *i nv（qn）; wn=fnn*inv（qn）;%E交多项式求根（适用于对称问题）%E交多项式求根（适用于非对称性问题）应用正交配置法求解以下等温球形催化剂颗粒内反应物浓度分布，其浓度分布的数学 模型为：（1）标准化令x r/R，y C/Cs代入微分方程及边界条件得：丄2 x2dy 36yx2 dx dxx 0単0dxx 1, y 1（2）离散化N 1Bji yi 36yj 0i 1j 1, 2,（

7、3）转化为代数方程组（以 N 3为例）r1 2 3 4 yyyy因为yN i目4 1，所以整理上式得:B11 36B12B13B21B22 36B23B31B32B33 36B41B42B43yi y2 y3B14B24B34B44 36本例中的代数方程组为线性方程组，可采用线性方程组的求解方法；若为非线性方程组则 采用相应的方法求解。N=3,权函数为1-xam,bm,wm,a n,bn,wn =collmatrix（3,3,1,3,1）; b1=bm;b1（i,i）=bm（i,i）-36;a0=b1（1:4,1:3）;b0=-b1（1:4,4）;y=aObO;y（4）=1;p=exam31（

8、3,3）;（注意要对文件修改权函数为 x=0.3631,0.6772,0.8998,1; % 零点ox0=0:0.1:1; % 真实解y0=si nh（6*x0）./x0/si nh（6）;（只用对称性配置矩阵）1-x2）若权函数改为1,则以下语句修改，其他不变am,bm,wm,a n,b n,w n=collmatrix（3,3,0,3,1）;（注意要对文件修改权函数为 1）x=0.4058,0.7415,0.9491,1;计算结果:权函数为1%0.9正交配置法0.80.70.6-y 0.50.40.3X：0.20.1ii1 0.50.7 0.8 0.9x边值问题的MatLab解法精确解:2

9、x（collfuni.m ） fun ctio n f=collfu ni（x,y） f（i）=y（2）;f（2）=4*y（1）; f=f（1）;f（2）;（collbci.m ） function f=collbc1（a,b） f=a（1）-1;b（1）-exp（2）;soli ni t=bvpi ni t（0:1,1,1） sol=bvp4c（collfu n1,collbc1,soli nit） plot（sol.x,sol.y）plot（sol.x,exp（2*sol.x）,*yi y2 真实解y x 1 y 2y 1 x2 e x 0 x 1y 0 2, y 1 1/eyi y2y2 1 x2 e x 2y1% 0 2, y1 1x 1 y21/e（collfun2.m ）function f=collfun2（x,y）f（2）=（1-x.A2）.*exp（-x）+2*y（1）-（x+1）.*y（2）; f=f（1）；（collbc2.m ）function f=collbc2（a,b） f=a（2）-2;b（2）-exp（-1）;soli nit=bvpi ni t（0:1,1,1）; sol=bvp4c（co