实验16常微分方程初值问题的数值解法 参考答案 1.利用Euler方法求初值问题的程序：function E=eulers(f,ab,yah)(见代码文件) 取步长分别为h=0.5和h=0.25求得的序列分别为 00.5115T, 25 y10.750.68750.7656250.94921875121191406251.533935546875 0 0.250.5 075 1 125 15 1 0.875 0.796875 0.759765620.758544920.7887268060.846385955 1875 640625 810547 1.75 2 2.25 2.5 2.753 0.9280877111.030826741.15197340 1.289226721.440573381.604251714 y334229 741745 399027 849149 74300500129 易见，该问题的精确解为)=3eP-2+,故3)=1.66939048044529。所以1=3处的最终全月 误差为E0(3),0.5)=0.135454933570289,E0(3),0.25=0.0651387664439962。得到的步进结 果图如下。可见，对同一种方法，采用较小的步长误差比较小。 05 25 2.改进的Euler方法求初值问题的程序：function IME=-impeuler(E,ab,yah)(见代码文件) 取步长分别为h=0.5和h=025求得的序 分别为 100.5 1 1.5 2 2.5 3 0.84370.8310546870.9305114746093 1.117587089538 1.373114913702 1.68212102632 75 57 01 97 0 0.25 0.5 0.75 1.25 1.5 0.83807373 0.814080715 0.8221962560.8586576320.920143066 0.8984375 046875 179443 369352 576069 258561 175 225 2.5 2.75 3 1.00372005 1.10679973 1.22709663 1.362593126 1.511507994 1.672268776 068139 224216 8%2003 28175 29561 21409 1=3处的最终全局误差为E0(3),0.5)=-0.0127305458844067
实验 16 常微分方程初值问题的数值解法 参考答案 1. 利用 Euler 方法求初值问题的程序：function E=eulers(f,a,b,ya,h)（见代码文件） 取步长分别为 h = 0.5 和 h = 0.25 求得的序列分别为 t 0 0.5 1 1.5 2 2.5 3 y 1 0.75 0.6875 0.765625 0.94921875 1.2119140625 1.533935546875 t 0 0.25 0.5 0.75 1 1.25 1.5 y 1 0.875 0.796875 0.75976562 5 0.75854492 1875 0.788726806 640625 0.846385955 810547 t 1.75 2 2.25 2.5 2.75 3 y 0.928087711 334229 1.03082674 741745 1.15197340 399027 1.28922672 849149 1.44057338 743005 1.604251714 00129 易见，该问题的精确解为 y=3e-t/2 -2+t，故 y(3)=1.66939048044529。所以 t =3 处的最终全局 误差为 E(y(3), 0.5)= 0.135454933570289, E(y(3), 0.25)= 0.0651387664439962。得到的步进结 果图如下。可见，对同一种方法，采用较小的步长误差比较小。 2. 改进的 Euler 方法求初值问题的程序：function IME=impeuler(f,a,b,ya,h)（见代码文件） 取步长分别为 h = 0.5 和 h = 0.25 求得的序列分别为 t 0 0.5 1 1.5 2 2.5 3 y 1 0.8437 5 0.831054687 5 0.9305114746093 75 1.117587089538 57 1.373114913702 01 1.68212102632 97 t 0 0.25 0.5 0.75 1 1.25 1.5 y 1 0.8984375 0.83807373 046875 0.814080715 179443 0.822196256 369352 0.858657632 576069 0.920143066 258561 t 1.75 2 2.25 2.5 2.75 3 y 1.00372005 068139 1.10679973 224216 1.22709663 862003 1.362593126 28175 1.511507994 29561 1.672268776 21409 t =3 处的最终全局误差为 E(y(3), 0.5)= -0.0127305458844067
E0(3),0.25)=-0.00287829576879939。得到的步进结果图如下。可见，改进的Euler方法比 原始Eul▣方法误差明显变小，原因是改进的Eul:方法的最终全局误差阶为O),而原始 的Euler方法的误差阶为OP)。 3.经典RK4方法求初值问题的程序：function RK4=rk4m(E,ab,ya,h)(见代码文件) 取步长分别为h=0.5和h=0.25求得的序列分别为 100.5 1 15 25 3 08364237808196247709 091714229539 11036825982 13595574922 16694307617 125 6558 5024 202 662 9916 0 025 05 075 1 125 15 .897491455 0.8364036680.811869580819594033 0.8557865510.917102058 1 078125 237229 242375 65079 933871 308097 1.75 2 2.25 2.5 2.75 3 100058853 1.103640815 1.223959876 1.359516811.508521142 1.669392747 y011477 76586 40518 679056 64965 88701 1=3处的最终全局误差为E03),0.5)4.02813538697977×10 E0(3),0.25)-226744172548976x106。得到的步进结果图如下。可见，经典4阶Runge--Luta 方法比Eulr方法误差明显变小，原因是RK4方法的最终全局误差阶为O的)，比上两题的 Eur方法的误差阶更高，但RK4方法每一步需要求4个函数值。 4.使用Malab内置的odc23和odc45函数计算得到的步进结果图如下。对于odc23函数 整个区间上的最大误差出现在=2.26处，误差为E=0.000163137202267238:对于ode45
E(y(3), 0.25)= -0.00287829576879939。得到的步进结果图如下。可见，改进的 Euler 方法比 原始 Euler 方法误差明显变小，原因是改进的 Euler 方法的最终全局误差阶为 O(h 3 )，而原始 的 Euler 方法的误差阶为 O(h 2 )。 3．经典 RK4 方法求初值问题的程序：function RK4=rk4m(f,a,b,ya,h)（见代码文件） 取步长分别为 h = 0.5 和 h = 0.25 求得的序列分别为 t 0 0.5 1 1.5 2 2.5 3 y 1 0.83642578 125 0.81962847709 6558 0.91714229539 5024 1.1036825982 2025 1.3595574922 6626 1.6694307617 9916 t 0 0.25 0.5 0.75 1 1.25 1.5 y 1 0.897491455 078125 0.836403668 237229 0.81186958 242375 0.819594033 650794 0.855786551 933871 0.917102058 308097 t 1.75 2 2.25 2.5 2.75 3 y 1.00058853 011477 1.103640815 76586 1.223959876 40518 1.35951681 679056 1.508521142 64965 1.669392747 88701 t =3 处的最终全局误差为 E(y(3), 0.5)= -4.02813538697977×10-5， E(y(3), 0.25)= -2.26744172548976×10-6。得到的步进结果图如下。可见，经典 4 阶 Runge-Lutta 方法比 Euler 方法误差明显变小，原因是 RK4 方法的最终全局误差阶为 O(h 4 )，比上两题的 Euler 方法的误差阶更高，但 RK4 方法每一步需要求 4 个函数值。 4．使用 Matlab 内置的 ode23 和 ode45 函数计算得到的步进结果图如下。对于 ode23 函数， 整个区间上的最大误差出现在 t=2.26 处，误差为 Err=0.000163137202267238；对于 ode45
函数，整个区间上的最大误差出现在0.15处，误差为Em=1.35671083256739×107。实际 上，这两个函数是不同阶的自适应步长的Runge--K方法：0c23使用的是显式的2阶-3 阶Rg-Ka方法：而od45则使用了显式的4阶-5阶ReeK方法。在每二步中 使用两个不同的求近似解的方法，并比较计算结果：如果两个结果的差小于给定误差限，则 接受该近似：如果两个结果的差大于给定误差限，则减小步长：如果答案超过了要求的有效 数字位数，则增加步长
函数，整个区间上的最大误差出现在 t=0.15 处，误差为 Err=1.35671083256739×10-7。实际 上，这两个函数是不同阶的自适应步长的 Runge-Kutta 方法：ode23 使用的是显式的 2 阶-3 阶 Runge-Kutta 方法；而 ode45 则使用了显式的 4 阶-5 阶 Runge-Kutta 方法。在每一步中， 使用两个不同的求近似解的方法，并比较计算结果：如果两个结果的差小于给定误差限，则 接受该近似；如果两个结果的差大于给定误差限，则减小步长；如果答案超过了要求的有效 数字位数，则增加步长
