常微分方程组(边值)

《常微分方程组(边值)》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、.-常微分方程组边值问题解法打靶法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(no

2、rm(c1-xn)>=eps&norm(x2-x1)>=eps)x0=x1;x1=x2;[a,b]=ode45(fun,[0,10],[0,x0]');c0=b(length(b),1);-.word.zl.-[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;应用打靶法求解以下边值问题：解：将其转化为常微分方程组的初值问题命令：x0=[0:0.1:10];-.word.zl.-y0=32*((cos(5)-1

3、)/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))-.word.zl.-函数：〔od

4、ebvp.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)-.word.zl.-计算结果：〔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’)-.word.zl.-有限差分法FiniteD

5、ifferenceMethodsFDM〔difference.m〕同上例：假设划分为10个区间，那么：-.word.zl.-函数：〔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=ab;-.word.zl.-x(

6、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,

7、y,’--’)-.word.zl.-[x,y]=difference(0,10,0,0,50);plot(x,y,’-.’)-.word.zl.-正交配置法OrthogonalCollocatioinMethodsCM构造正交矩阵函数〔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

8、=1:mxm(i)=x0(m+1-i);endxm(m+1)=1;forj=1:m+1fori=1:m+1-.word.zl.-qm(j,i)=xm(j)^(2*i