常微分组边值问题打靶法matlab解法.docx

常微分方程组两点边值问题的数值解法----张亚苗2011年9月y 〃 —4y = 0 勺=：2y = 4 y y(0) = 1 可化为微分方程组1 y()=1 y(1) = 3y⑴=3方法一：配置法Matlab 程序： function bvcollation clcsolinit = bvpinit(linspace(0,1,20),[100 600]);% sol = bvp4c(@twoode,@twobc,solinit);x = linspace(0,1,20);y = deval(sol,x);y'plot(x,y(1,:),x,y(2,:));end%微分方程组function dydx = twoode(x,y)dydx = [ y(2)4*y(1)];end%边值条件function res = twobc(ya,yb)res = [ ya(1)-1yb(1)-3];end运行结果：1.0000 -0.42030.9834 -0.21170.9777 -0.00550.9828 0.20070.9988 0.40911.0259 0.62201.0644 0.84191.1147 1.07101.1774 1.31211.2531 1.56771.3427 1.84071.4472 2.13411.5678 2.45121.7057 2.79541.8626 3.17072.0401 3.58112.2402 4.03132.4652 4.52612.7175 5.07123.0000 5.67246方法二：打靶法程序：积分用ode45，搜索用二分法 function shoot_first_try clcclear alla=0;b=1;a2=-1;b2 = 0;s = 0;tol=1.0*10e-6;h1=Fun(a2)h2=Fun(b2)s=bisect(@Fun,a2,b2,tol)[t1,y1]=ode45(@f,[a,b],[1,s]) plot(t1,y1(:,1),t1,y1(:,2)) function xc=bisect(f,a,b,tol)if sign(f(a))*sign(f(b))>=0 error('f(a)*f(b)<0 not satisfied!') end fa=f(a);fb=f(b);k=0;while (b-a)/2>tol c=(a+b)/2;fc=f(c);if fc==0breakendif sign(f(a))*sign(f(b))<0 b=c;fb=fc;elsea=c;fa=fc;endendxc=(a+b)/2;function z=Fun(s)a=0;b=1;yb=3;h=linspace(a,b,20);y0=[1;s];[t,y]=ode45(@f,h,y0);z=y(end,1)-yb;function ydot=f(t,y)%ydot=[0;0];ydot(1)=y(2);ydot(2)=4火y(1);ydot=[ydot(1);ydot(2)];运行结果：1.0000 -0.50000.9969 -0.47490.9940 -0.44990.9913 -0.42500.9887 -0.40010.9799 -0.30170.9736 -0.20410.9697 -0.10690.9683 -0.01000.9692 0.08680.9726 0.18390.9784 0.28140.9867 0.37970.9974 0.47881.0106 0.57921.0264 0.68111.0447 0.78461.0656 0.89011.0892 0.99781.1155 1.10801.1446 1.22101.1766 1.33701.2115 1.45641.2495 1.57941.2905 1.70641.3348 1.83771.3825 1.97351.4335 2.11431.4882 2.26031.5466 2.41201.6089 2.56981.6751 2.73391.7456 2.90491.8204 3.08321.8998 3.26921.9840 3.46332.0731 3.66612.1674 3.87812.2671 4.09982.3724 4.33172.4837 4.57452.5711 4.76372.6621 4.95962.7569 5.16252.8555 5.3726。