常微分組邊值問題打靶法matlab解法

1、常微分方程組兩點邊值問題的數(shù)值解法 -張亞苗2011年9月可化為微分方程組方法一：配置法Matlab程序：function bvcollationclc solinit = 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 re

2、s = twobc(ya,yb) res = ya(1)-1 yb(1)-3;end運行結(jié)果： 1.0000 -0.4203 0.9834 -0.2117 0.9777 -0.0055 0.9828 0.2007 0.9988 0.4091 1.0259 0.6220 1.0644 0.8419 1.1147 1.0710 1.1774 1.3121 1.2531 1.5677 1.3427 1.8407 1.4472 2.1341 1.5678 2.4512 1.7057 2.7954 1.8626 3.1707 2.0401 3.5811 2.2402 4.0313 2.4652 4.52

3、61 2.7175 5.07123.0000 5.6724方法二：打靶法程序：積分用ode45,搜索用二分法function shoot_first_tryclcclear 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)<

4、;0 not satisfied!')endfa=f(a);fb=f(b);k=0;while (b-a)/2>tol c=(a+b)/2; fc=f(c); if fc=0 break end if sign(f(a)*sign(f(b)<0 b=c;fb=fc; else a=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;ydo

5、t(1)=y(2);ydot(2)=4*y(1);ydot=ydot(1);ydot(2); 運行結(jié)果： 1.0000 -0.5000 0.9969 -0.4749 0.9940 -0.4499 0.9913 -0.4250 0.9887 -0.4001 0.9799 -0.3017 0.9736 -0.2041 0.9697 -0.1069 0.9683 -0.0100 0.9692 0.0868 0.9726 0.1839 0.9784 0.2814 0.9867 0.3797 0.9974 0.4788 1.0106 0.5792 1.0264 0.6811 1.0447 0.7846 1.0656 0.8901 1.0892 0.9978 1.1155 1.1080 1.1446 1.2210 1.1766 1.3370 1.2115 1.4564 1.2495 1.5794 1.2905 1.7064 1.3348 1.8377 1.3825 1.9735 1.4335 2.1143 1.4882 2.2603 1.5466 2.4120 1.6089 2.5698 1.6751 2.7339 1.7456 2.9049 1.8204 3.0832 1.8998 3.2692 1.9840 3.4633 2.0731 3.6661 2.1674 3.8781 2.2