1、数学基础：微分方程（组）数值解法基于的实现与分析1 常微分方程初值问题的数值解法微分方程（组）是科学研究和工程应用中最常用的数学模型之一。如揭示质点运动规律的 Newton 第二定律：（）002xtxtFdm和刻画回路电流或电压变化规律的基尔霍夫回路定律等，但是，只有一些简单的和特殊的常微分方程及常微分方程组，可以求得用公式给出的所谓“解析解”或“公式解”，如一阶线性微分方程的初值问题：（）0ytfadt的解为：（）dfetytaat0但是，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，因此，基于如下的事实：、绝大多数的常微分方程和常微分方程组得不到解析解；、实际应用中往往只
2、需要知道常微分方程（组）的解在（人们所关心的）某些点处的函数值（可以是满足一定精度要求的近似值）；如果只需要常微分方程（组）的解在某些点处的函数值，则没有必要非得通过求得公式解，然后再计算出函数值不可，事实上，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。一般的一阶常微分方程初值问题是指如下的一阶常微分方程的定解问题：（）0,ytttfdf微分方程（组）的初值问题通常是对一动态过程演化规律的描述，求解常微分方程（组）的初值问题就是要了解和掌握动态过程演化规律。一、最简单的数值解法Euler 方法假设要求在点（时刻） ， ， 处初值问khtk0ntf0nk,2
3、1题（）的解 的近似值。首先对式（）的两端积分，得ty（）dtyftydt kk tkt 11 ,1对于（）的右边，如果被积函数用积分下限 处的函数值代替被k积函数作积分（从几何上的角度看，是用矩形面积代替曲边梯形面积），则有（）ktkk tyhfdtyfytk ,11 进而得到下式给出的递推算法Euler 方法（）01 ,21,yt nkthfkk例 用 Euler 方法解如下初值问题，取 ，3.0h10sin23yttdt解：由（）得1010,2sin6.3.3yktkkk 结果如下：Euler_Method 0 5 1 5 2 5 37 8 9 11 2 4 m 如果取 ，其结果如下图所
4、示：.hEuler_Method 5 5 5 37 8 11 2 3 4 附 应用程序 Euler_Method.m%Euler Method for the Example: y = y - 2sin(t)y3; y(0)=1;clear allT,z=ode45(Euler_Method_file,0;3,1);t=0:0.3:3;%t=0:0.1:3;%t=0:0.3:3;Lt=length(t);y=zeros(1,Lt);y(1)=1;for k=2:Lty(k)=1.3*y(k-1)-0.6*sin(t(k-1)*y(k-1)3;%y(k)=1.1*y(k-1)-0.2*sin(t
5、(k-1)*y(k-1)3;%y(k)=1.3*y(k-1)-0.6*sin(t(k-1)*y(k-1)3;endplot(T,z,r)hold onplot(t,y,t,y,r*)title(Euler Method of Solving Initial Value Problem)legend(Integral Curve,Euler Curve)附 应用程序 Euler_Method_File.mfunction varargout = odefile1(t,y,flag)switch flagcase % Return dy/dt = f(t,y).varargout1 = f(t,y
6、);case init % Return default tspan,y0,options.varargout1:3 = init;case jacobian % Return Jacobian matrix df/dy.varargout1 = jacobian(t,y);otherwiseerror(Unknown flag flag .);end% -function dydt = f(t,y)dydt =y(:)-2.*sin(t(:).*y(:).3;% -function tspan,y0,options = inittspan =0;3; y0 = 1;options = ode
7、set(jacobian,on);% -%function dfdy = jacobian(t,y)dfdy =1-2.*(cos(y(:).*y(:).3+3.*sin(t(:).*y(:).2);% -二、改进的 Euler 方法对于（）的右边，如果被积函数用积分限 和 处的函数值的kt1kt算术平均值代替（几何上，是用梯形面积代替曲边梯形面积），则有（）11 ,2,1 kktkk tyftyfhdtyfytk进而得到下式给出的递推算法：（）0 11 ,2,2yt nytfytfhkkk 通常算法（）比 Euler 方法（）的精度高，但是，按算法（）求 时要解（非线性）方程，这是算法（）不
8、如 Euler 1k方法的方面，为了） 尽可能地保持算法（）精度高的优点；） 尽可能地利用 Euler 方法计算简单的长处；人们采取了如下的称之为改进的 Euler 方法的折衷方案：预测 kkythfy,01修正 （）nkytfkkk ,21,2010yt例 Euler 方法与改进的 Euler 方法的比较下图是当 时比较的结果：3.hImproved_Euler_Method 0 5 1 5 2 5 3 8 1 2 3 4 d 附 应用程序 Improved_Euler_Method.m%Compare Euler Method With Improved Euler Methodclear
9、 allpackclfT,z=ode45(Euler_Method_file,0;3,1);t=0:0.3:3;Lt=length(t);y=zeros(1,Lt);y(1)=1;x=zeros(1,Lt);x(1)=1;for k=2:Lty(k)=1.3*y(k-1)-0.6*sin(t(k-1)*y(k-1)3;x(k)=x(k-1)+0.15*(x(k-1)-2*sin(t(k-1)*x(k-1)3+y(k)-2*sin(t(k)*y(k)3);endplot(T,z,r)hold onplot(t,y,t,x,g)plot(t,y,r*,t,x,*)title(Compare Eul
10、er Method With Improved Euler Method)legend(Integral Curve,Euler Curve,improved Euler Curve)三、Euler 方法和改进的 Euler 方法的误差分析由 Taylor 公式（14）22111, !hOtyfty ttytkk kkk 说明 Euler 方法的截断误差是 ，类似地，由（15）321 !1, hOtytftyt kkkk （16）3211,thtftt kkkk 以及Otytkk1可得（17）311 ,2hOtyftyfhtyt kkkk 从上述推导 Euler 方法、改进的 Euler 方法
11、的过程以及例、例容易看出，改进的 Euler 方法 Euler 方法的精度高，其原因在于：） 在推导 Euler 方法时，我们是用待求解函数 在一点处的ty变化率 代替 在区间 上的平均变化率；ktyf,ty,1kt） 而在推导改进的 Euler 方法时，我们是用待求解函数 在t两点处变化率的平均值 代替 在区间1,21kktyftyf ty上的平均变化率；,1kt显然，通常 比 更接近于 在1,21kktyftyf ktyf, ty区间 上的平均变化率。由此，人们受到启发：适当地选取区间,1kt上函数 多点处的变化率，用它们加权平均值代替 在ktty ty区间 上的平均变化率，近似解的精度应
12、更高。,1kt下面将要介绍的 RungeKutta 法就是基于上述想法得到的。四、RungeKutta 法RungeKutta 法是按选取区间 上函数 变化率 的个,1kttyytf,数的多少和截断误差的阶数 来区分的一系列方法，如mhO） 二阶的 RungeKutta 法（改进的 Euler 方法）（18）21121,KhytfKkk） 三阶的 RungeKutta 法（19）3211312146,KhytfKytfkkk） 四阶的 RungeKutta 法（） 古典形式（20）432114231216,KKhytfhytfKkkkk（） Gill 公式(具有减小舍入误差的优点)（21） 432114 213216,2, KKhyhtfKyhtfKtfKkkkkkk例三：疾病的传播与防疫模型其中 表示一种传染病传染者的人数， 表示易感染者人数。tx tyedxycaxdtyb2Runge_Kutta_Method 0 8 5 例 2：二体问题Runge_Kutta_Method01
