龙格席塔方法的 MATLAB实親 实例1海上私() ()=f(t,x),x(t0)=x0,x=(x1,…xn),f=(f1,…fn)y 模型的数值解 [x=0de23(@t,ts,xo)32阶龙格库塔公式 x=ode45(@t,ts,x)5氟4阶龙格库塔公式 +(ar 是待解方程写成 function dx=f(tx) xin,(o) 的函数m文件 √(c-x)2+(at-y ts=[to,t1,…,t输出指定时刻to,t,….,t的函数值 ts to: k: tf 输出[to切内等分点处的函数值 船速a=20(海里小时) 求:私晨的世置x(2y(0 艇速b=40(海里/小时) 私的航能y(x) X0为函数初值(n维)输出t=ts,x为相应函数值(n维) 缺省精度(相对误差103,绝对误差106, 距高c=15(海里 计算步长按精度要求自动调整 output t,x(D), y(t) ∥寥例1海上私()型的激 a=20b=40,c=15; a=20,b=40,c=15 x =ode 45(@jisi, ts, xo, D, a, b, c) plot(t, x),grid, gtext(,x(t), FontSize, 16 走私船的位量 gtext('y(ty, FontSize, 16),pause x1(D)=c=15 0103.98540.29242.0 %draw y(x): the position of match js %Creat the function for jisi. m plot(x( 1), x( 2), T"), grid 055945410690630 缉私艇的航线v(x) 0207851519940 sqrt(c-x(1)^2+a·-x(2)2 0.301134963.20056.0 dx=|b°(c-x(1)/sb°a·tx(2)s 0351217045570 =0.5时缉私艇追上走私船0451474518.027390 被的值解 学学实 实创1 模型的数值解学醇学实妇 实1 海上() yoo v,a 海上輯私() 设,不变,凌变大03%58 为30,35接近40,03105240482831105 观察解的变化 041253482755140 地的航yx 0412539454826984014 a=35,b=4,c=150513751120830175 0.51375397412.075344175 ?缛私艇追上走私船□121498940144:20 214996100010 13「14996441514551 13|149941001455 积慑戆大 1415011480183490 提高精度 产6时缉私艇追上走私船141900 515003520146525 6149865594865 判“地上的有欺亦海?|Ld1sms9g3 龙格—库塔方法的 MATLAB 实现 T n T n x(t) f (t, x), x(t ) x , x (x , x ) , f ( f , f ) & = 0 = 0 = 1 L = 1 L [t,x]=ode23(@f, ts, x0) 3级2阶龙格-库塔公式 [t,x]=ode45(@f, ts, x0) 5级4阶龙格-库塔公式 f是待解方程写成 的函数m文件： function dx=f(t,x) dx=[f1; f2;…; fn]; ts = [t0,t1, …,tf] 输出指定时刻 t0,t1, …,tf 的函数值 ts = t0:k:tf 输出 [t0,tf] 内等分点处的函数值 x0为函数初值(n维) 输出t=ts, x为相应函数值(n维) 缺省精度(相对误差10-3，绝对误差10-6), 计算步长按精度要求自动调整. 实例1 海上缉私(续) 模型的数值解 2 2 2 2 ( ) ( )( ) ( ) ( )( ) dx b c x dt c x at y dy b at y dt c x at y ⎧ − = ⎪ ⎪ −+− ⎨ − ⎪ = ⎪ − +− ⎩ x y (0) 0, (0) 0 = = 0 y x c (x(t),y(t)) a b 设: 船速a=20 (海里/小时) 艇速b=40 (海里/小时) 距离c=15 (海里) 求: 缉私艇的位置x(t),y(t) 缉私艇的航线y(x) x0=[0 0]; a=20;b=40;c=15; [t,x]=ode45(@jisi,ts,x0,[],a,b,c); %exact solution x1=c y1=a*t; %output t,x(t),y(t) [t,x,y1] %draw x(t),y(t) plot(t,x),grid, gtext('x(t)','FontSize',16), gtext('y(t)','FontSize',16),pause %draw y(x): the position of tatch js plot(x(:,1),x(:,2),'r*'),grid xlabel('x','FontSize',16), ylabel('y','FontSize',16) %Creat the function for jisi.m %Let x(1)=x, x(2)=y function dx=jisi(t,x,a,b,c) s=sqrt((c-x(1))^2+(a*t-x(2))^2); dx=[b*(c-x(1))/s;b*(a*t-x(2))/s]; MATLAB 6.5.1.lnk jisi.m, seajisi.m 实例1 海上缉私(续) 模型的数值解 a=20, b=40, c=15 走私船的位置 x1(t)= c=15 y1(t)=at=20t t=0.5时缉私艇追上走私船 0 5 10 15 20 0 2 4 6 8 10 x y 缉私艇的航线y(x) 0.50 15.0046 9.9979 0.45 14.7451 8.0273 0.40 13.9806 6.1773 0.35 12.8170 4.5552 0.30 11.3496 3.2005 0.25 9.6705 2.1178 0.20 7.8515 1.2899 0.15 5.9445 0.6906 0.10 3.9854 0.2924 0.05 1.9984 0.0698 0 0 0 t x(t) y(t) 10.0 9.0 8.0 7.0 6.0 5.0 4.0 3.0 2.0 1.0 0 y1(t) 实例1 海上缉私(续) 模型的数值解 设b，c不变，a变大 为30, 35, …接近40, 观察解的变化 : a=35, b=40, c=15 t=? 缉私艇追上走私船 1.6 14.9866 55.9486 56.0 1.5 15.0023 52.0146 52.5 1.4 15.0117 48.0183 49.0 1.3 14.9996 44.0165 45.5 1.2 14.9986 40.0164 42.0 … … … … 0.5 13.7551 12.0830 17.5 0.4 12.5384 8.2755 14.0 0.3 10.5240 4.8283 10.5 0.2 7.5928 2.1308 7.0 0.1 3.9561 0.5058 3.5 0 0 0 0 y1 t x(t) y(t) (t) 累 积误差较大 提高精度! 实例1 海上缉私(续) 模型的数值解 a=35, b=40, c=15 opt=odeset('RelTol',1e-6, 'AbsTol',1e-9); [t,x]=ode45(@jisi,ts,x0,opt); t=1.6时缉私艇追上走私船 0 5 10 15 20 0 20 40 60 x y 缉私艇的航线y(x) 判断“追上”的有效方法? 1.6 15.000020 55.999931 56.0 1.5 14.999998 52.000005 52.5 1.4 14.999993 48.000005 49.0 1.3 14.999963 44.000005 45.5 1.2 14.999616 40.000005 42.0 … … … … 0.5 13.753974 12.075344 17.5 0.4 12.539454 8.269840 14.0 0.3 10.521921 4.829308 10.5 0.2 7.592822 2.130678 7.0 0.1 3.956104 0.505813 3.5 0 0 0 0 y1 t x(t) y(t) (t)