常微分方程组(边值).doc

PAGE 常微分方程组边值问题解法 打靶法Shooting Method（shooting.m） %打靶法求常微分方程的边值问题 function [x,a,b,n]=shooting(fun,x0,xn,eps) if nargin3 eps=1e-3; end x1=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; end x=x2; 应用打靶法求解下列边值问题： 解：将其转化为常微分方程组的初值问题 命令： 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.9524 plot(x,y(:,1)) x0=[0:1:10]; y0=32*((cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); hold on plot(x0,y0,’o’) 有限差分法Finite Difference Methods FDM（difference.m） 同上例： 若划分为10个区间，则： 函数：（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-2 a(i,i+1)=1; a(i+1,i)=1; end b=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:n x(i)=x0+(i-1)*h; y(i)=yy(i-1); end x(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:m xm(i)=x0(