完整版常微分方程组边值.docVIP

. 常微分方程组边值问题解法 打靶法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; 应用打靶法求解下列边值问题： 2?yyd?8?? 24dx????00y ?????y010?? 解：将其转化为常微分方程组的初值问题 教育资料． . dy?1?y?2dx?dyy?12??8? dx4????0yy?01?dy?t?dx??0x 命令： 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） 同上例： 2y?yy?2yydyii1?1ii??8???8? 2244dxh2??h2??h8???2yy?y ??1i?ii1? 4?? 10个区间，则：若划分为??2??h??1?2????? 4??????2??hy????2????21?1?0?8h1?? 4??????y2??h82??????????? ???????2yh82??hn?2????????2??11208h?y?????? ??1n?4????2??h????2?1????? 4???? 函数：（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