081数值计算方法—常微分方程(组).doc

科学计算理论、方法及其基于MATLAB的程序实现与分析微分方程（组）数值解法1 常微分方程初值问题的数值解法微分方程（组）是科学研究和工程应用中最常用的数学模型之一。如揭示质点运动规律的Newton第二定律：（）和刻画回路电流或电压变化规律的基尔霍夫回路定律等，但是，只有一些简单的和特殊的常微分方程及常微分方程组，可以求得用公式给出的所谓“解析解”或“公式解”，如一阶线性微分方程的初值问题：（）的解为：（）但是，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，因此，基于如下的事实：、绝大多数的常微分方程和常微分方程组得不到（有限形式的）解析解；、实际应用中往往只需要知道常微分方程（组）的解在（人们所关心的）某些点处的函数值（可以是满足一定精度要求的近似值）；如果只需要常微分方程（组）的解在某些点处的函数值，则没有必要非得通过求得公式解，然后再计算出函数值不可，事实上，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。一般的一阶常微分方程（组）的初值问题是指如下的一阶常微分方程（组）的定解问题：（）其中（） （）常微分方程（组）的初值问题通常是对一动态过程（动态系统、动力系统）演化规律的描述，求解常微分方程（组）的初值问题就是要了解和掌握动态过程演化规律。.1 常微分方程（组）的Cauch问题数值解法概论假设要求在点（时刻），处初值问题（）的解的（近似）值,如果已求得时刻的值或它的近似值(如时刻的值)，那么将式（）的两端在区间上积分(10)可得 (11)或(12)显然,为了利用式（11）或（12）求得的精确值（近似值），必须计算右端的积分，这是问题的关键也是难点所在，如前所述，一般得不到精确的公式解，因此需要采用数值积分的方法求其近似解,可以说,不同的式值积分方法将给出不同的Cauch问题的数值解法。1.2最简单的数值解法Euler 方法假设要求在点（时刻），处初值问题（）的解的近似值。首先对式（）的两端积分，得（13）对于式（13）的右边，如果用积分下限处的函数值代替被积函数作积分（从几何上的角度看，是用矩形面积代替曲边梯形面积），则有（14）进而得到下式给出的递推算法Euler 方法（15）例用Euler 方法解如下初值问题，取，解：由（15）得结果如下：open Euler_Method.m 如果取，其结果如下图所示：Euler_Method 1.3改进的Euler 方法 对于（15）的右边，如果被积函数用积分限和处的函数值的算术平均值代替（几何上，是用梯形面积代替曲边梯形面积），则有（16）进而得到下式给出的递推算法：（17）通常算法（17）比Euler 方法（15）的精度高，但是，按算法（17）求时要解（非线性）方程（组），这是算法（17）不如Euler 方法的方面，为了） 尽可能地保持算法（17）精度高的优点；） 尽可能地利用Euler 方法计算简单的长处；人们采取了如下的称之为改进的Euler 方法的折衷方案：预测 （18）修正 （19） 例Euler 方法与改进的Euler 方法的比较下图是当时比较的结果：open Improved_Euler_Method.m 1.4Euler 方法和改进的Euler 方法的误差分析由Taylor 公式 （19）说明Euler 方法的截断误差是，类似地，由 （20） （21）以及 (22)让式（20）的两端减式（21）的两端，可得 （23）从上述推导Euler 方法、改进的Euler 方法的过程以及例、例容易看出，改进的Euler 方法Euler 方法的精度高，其原因在于：1 在推导Euler 方法时，我们是用待求解函数在一点处的变化率代替在区间上的平均变化率： （24）2 而在推导改进的Euler 方法时，我们是用待求解函数在两点处变化率的平均值代替在区间上的平均变化率；显然，通常比更接近于在区间上的平均变化率。由此启发人们：适当地选取区间上函数若干点处的变化率，用它们加权平均值代替在区间上的平均变化率，近似解的精度应更高。下面将要介绍的RungeKutta法就是基于上述想法得到的。2RungeKutta法RungeKutta法是按选取区间上函数变化率的个数的多少和截断误差的阶数来区分的一系列方法，如 二阶的RungeKutta法（改进的Euler 方法） （25） 三阶的RungeKutta法 （26） 四阶的RungeKutta法 1) 古典形式 （27） 2) Gill 公式(具有减小舍入误差的优点) （28）4 RungeKutta法的一般形式 (29)其中称为增量函数(Increment Function)以及 （30）需要特别指出的是：在确定（29）、（30）中的参数，和时，应该使（29）的右端和适当阶的（如阶）Taylor展式一致，这样，至少对于底阶的R-K法来说，参与加权的斜率个数与方法的阶数是一致的。例如：一阶的R-K法即Euler法，局部的截断误差（Truncation Error）是，所以，当微分方程的阶是一次函数时是精确的；二阶R-K法即改进的Euler法，局部的截断误差是，所以，当微分方程的阶是二次函数时是精确的；类似地，四阶R-K法局部的截断误差是。但是，一般五阶以及五阶以上的R-K法的局部截断误差不再具有上述的特点，因此用“性价比”来衡量，四阶R-K得到了广泛的应用。下面是Butchers Fifth-Order R-K 方法（1964）： （28） （29）例三：疾病的传播与防疫模型open Lunge_Kutta_Method.m其中表示一种传染病传染者的人数，表示易感染者人数。 (31)该系统的Jacobi矩阵为： (32)该系统的平衡点为： (33)根据问题的实际意义，有意义的平衡点为： (34)3解刚性问题的数值方法例四：单个微分方程的例子考虑如下的一阶微分方程描述的系统 (35)容易求的该微分方程初值问题的通解： (36)f=dsolve(Dy = -1000*y+3000-2000*exp(-t), y(0) = c) 并且有: (37)当初值时，方程的解中包含“快变”和“慢变”两部分.应用不同数值方法求解的结果：解非刚性问题的方法：ODE45: Sol_StiffODE01.m解刚性问题的方法：ODE15s: Sol_StiffODE01.m解非刚性问题的方法：Euler Method: Sol_StiffODE01_Euler.m系统（35）所对应的自治系统为 (38)利用Euler 法求解上述初值问题：， (39)为保证时间序列的收敛性，步长应满足不等式：.例五：微分方程组的例子考虑如下的一阶微分方程组描述的系统 (40) (41) (42)利用Euler 法求解上述初值问题的迭代公式为， (43)由于系数矩阵的特征值为(44)所以Euler 法不适用于该初值问题的求解。 利用四阶R-K 法求解上述初值问题的迭代公式为， (45)其中系数矩阵为(46)其特征值为： (47)解非刚性问题的方法：ODE45: Sol_StiffODE02.m解刚性问题的方法：ODE15s: Sol_StiffODE02.m解非刚性问题的方法：Euler Method: Sol_StiffODE02_Euler.m解非刚性问题的方法：Euler Method: Sol_StiffODE02_RK4.mf=dsolve(Dx=-5*x+3*y,Dy=100*x-301*y,x(0)=1,y(0)=2)transient dominate motivationorientationscopediscretizationpredictor corrector approachincrement function recurrenceadaptive steps - size controlabrupt change例2：二体问题open Runge_Kutta_Method014 线性多步法与常微分方程数值解法的分类 一般的常微分方程初值问题的数值解法都是以递推的形式给出的，即是递推算法，根据递推算法，在计算时，已经得到了前面各时刻的近似值，前面介绍的各种数值解法都有一个共同点：在计算时，只用到了前一个时刻（当前时刻）的“信息”预测未来某时刻系统的“状态”，这样的数值解法称为单步法，对于一个动态过程在时刻的状态而言，不仅前一个时刻的信息对它有影响，前若干个时刻的信息通常对它也有影响，显然，单步法的缺点是没有充分利用已得到的信息。基于上述考虑，适当取前若干个时刻的信息，并用的线性组合代替在区间上的平均变化率所给出的方法称为线性多步法。4.1 线性多步法: 从二步法谈起 在改进的Euler方法中预测 （48）较正 （49） 如果改进预测的手段，如利用数值微分的研究成果： (50)那么，由于（48）的局部截断误差是，精度高于Euler公式的局部截断误差，所以，下列预测-校正方法Predictor: （51）Corrector: （52） 尽管步长增大了，精度却提高了。公式（49）、（50）给出的数值方法称为非自始的（non-self-starting）二步法，其原因在于：利用该方法求解微分方程初值问题时，自身不能提供所需的初值。4.3 线性二步法的进一步改进在上述讨论中，我们只考虑了预测手段方面的改进，从方法论的角度进行思考，也应该考虑校正方面的改进问题，并且，可以希望适当的改进会带来精度提高的回报。下面推荐的方案就是常用的一种：Predictor: （53）Corrector: , （54）其中是时刻最终计算结果,校正的次数按如下公式确定 （55）其中是预先指定的误差精度（Prespecified Error Tolarence）.4.4 线性二步法的误差分析与步长控制1 预测误差： (56)其中下表表示预测误差，这表明 （57）2 校正误差：类似地可得 （58）3 步长控制准则有误差估计（53）-（55）得到 （59）如果大于事先指定的误差精度，那么就减小步长。4.5 线性多步法的进一步扩展 在积分公式 （60）中，用个节点，的次插值多项式代替被积函数，得到的就是一般的线性多步法。如当时， （61）如当时， （62）特别地，用四个节点，的三次插值多项式代替被积函数，得到 （63）四个节点，的三次插值多项式代替，又得到 （64）其中误差分别为、，可见，公式（64）的精度比公式（63）的精度高。Adams 预报-校正公式： （65）其中，说明：公式（65）自身不能提供最初所需的三个初值，需要Runge-Kutta 法的配合。5 应用常微分方程组初值问题数值解法的基本要求5.1 常微分方程组初值问题数值解法的本质以Runge-Kutta为例，对于给定的常微分方程（组）的初值问题：（）其一般形式是 (29) 从一般性的系统角度看，（7）刻画的是连续时间动态（力）系统的演化规律，而（29）刻画的是离散时间动态系统的演化规律，由此可见，常微分方程（组）初值问题数值解法的本质是用离散时间动态系统模型代替连续时间动态系统模型，这意味着，应用数值方法求解对于给定的常微分方程（组）的初值问题时，所选择的离散时间动态系统模型应该是连续时间动态系统模型的“有意义的近似”。 为此，应用常微分方程组初值问题数值解法需要了解的三个基本概念：1 相容性：因为 (66)所以，应该保证 (67)2 收敛性：如果已知时刻函数的精确值， (68)那么应该有 (69)3 稳定性：5 常微分方程边界值问题的数值解法对于给定的阶常微分方程组(70)其定解问题要求给定个独立的定解条件,如果个独立的定解条件指定在求解区间的一端,那么这样的定解问题就是前面讨论过的初值问题,如果个独立的定解条件分别在求解区间的两端给定,这样的定解问题就是所谓的边界值问题(Two-Points Boundary Valve Problems).常微分方程的边界值问题在工程实际中广泛存在，如两端固定的弹簧振动问题，与外界有热交换的杆状物体的热平衡问题等。下面介绍两种常用的常微分方程边界值问题的数值解法。5.1 边界值问题的数值解法一: Shooting Method 采用 Shooting Method解常微分方程边界值问题的主导思想是将边界值问题转化为等价的初值问题。例6 考虑如下的边界值问题： （71）其中是细长杆状物体的长度，时周围空气的温度，时热交换系数（单位：），那么（71）就是描述该物体处于稳态时温度变化规律的数学模型，如果在杆的两端，外界的温度固定： （72）那么式（71）和（72）就构成了简单的常微分方程边界值问题，设参数和边界条件如下： （73）下面讨论该边界值问题数值解，首先将方程（71）改写成等价的一阶微分方程组的形式，令， （74）如果我们把问题（74）作为初值问题求解，那么缺少未知函数的初始条件，为此，先对它的值作出猜测，例如，然后选择适当的数值方法（如四阶Runge-Kutta法）进行求解，结果如下open Sol_BVProblem01.m根据上述计算结果，减小的猜测值为，计算结果为由于（74）是线性微分方程组，其公式解为 （75）（75）表明与之间是线性关系，从而可利用下述公式计算5.2 边界值问题的数值解法二: Finite-Differance Methods对于非线性边界值问题,由于初值与终值之间的关系是非线性，所以采用Shooting Method求解会遇到困难,且不易得到精确的结果，在这种情况下，可考虑采用有限差分方法续例6 利用数值微分的结果：二阶中心有限差商公式: (76)得到 (77)当取时，有（77）可得 （78）A=-2.04 1 0 0; 1 -2.04 1 0; 0 1 -2.04 1; 0 0 1 -2.04;b=-40.8; -0.8; -0.8; -200.8;T=Ab T = 65.9698 93.7785 124.5382 159.4795 附应用程序1: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:Lt y(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(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);case init % Return default tspan,y0,options. varargout1:3 = init;case jacobian % Return Jacobian matrix df/dy. varargout1 = jacobian(t,y);otherwise error(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 = odeset(jacobian,on);% -%function dfdy = jacobian(t,y)dfdy =1-2.*(cos(y(:).*y(:).3+3.*sin(t(:).*y(:).2);% -附应用程序Improved_Euler_Method.m%Compare Euler Method With Improved Euler Methodclear 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:Lt y(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 Euler Method With Improved Euler Method)legend(Integral Curve,Euler Curve,improved Euler Curve)附应用程序Runge_Kutta_Method.m%4th Order Runge_Kutta Method For Systems Equations%dx/dt=axy-b*x%dy/dt=-axy-cx2y+dx+eclear%,a,b,c,d are belong to interval (0,1)a=0.5;b=0.9;c=0.9;d=0.1;e=1;flag=jacobian;T,z=ode15s(Runge_Kutta,0;30,1;3,a,b,c,d,e);plot(T,z)%4th Order Runge_Kutta Method For Systems of Differential Equationsh=0.5;t=0:h:30;y=;z=;K1=;K2=;K3=;K4=;y(:,1)=1;3;for k=1:length(t)-1; K1(:,k)=a*y(1,k).*y(2,k)-b*y(1,k) -a*y(1,k).*y(2,k)-c*(y(1,k).2).*y(2,k)+d*y(1,k)+e; z(:,k)=y(:,k)+2h.*K1(:,k); K2(:,k)=a*z(1,k).*z(2,k)-b*z(1,k) -a*z(1,k).*z(2,k)-c*(z(1,k).2).*z(2,k)+d*z(1,k)+e; z(:,k)=y(:,k)+2h.*K2(:,k); K3(:,k)=a*z(1,k).*z(2,k)-b*z(1,k) -a*z(1,k).*z(2,k)-c*(z(1,k).2).*z(2,k)+d*z(1,k)+e; z(:,k)=y(:,k)+h.*K3(:,k); K4(:,k)=a*z(1,k).*z(2,k)-b*z(1,k) -a*z(1,k).*z(2,k)-c*(z(1,k).2).*z(2,k)+d*z(1,k)+e; y(:,k+1)=y(:,k)+8h*(K1(:,k)+3*(K2(:,k)+K3(:,k)+K4(:,k);endhold onplot(t,y(1,:),r*,t,y(2,:),g*)title(系统演化曲线)legend(Numbers of Patient,Numbers of Easy to be Infected,RK-Numbers of Patient,. RK-Numbers of Easy to be Infected,1)syms x y u v %a b c d e;x,y=solve(a*x*y-b*x=0,-a*x*y-c*x2*y+d*x+e=0); % find equilibrium J=jacobian(a*u*v-b*u;-a*u*v-c*u2*v+d*u+e,u v); % Evaluate JACOBI MatrixEJ=eig(J); % Evaluate Eigenvalues of Jacobi Matrixpausehold offequilibrium =1/2/c/b*(-b*a+d*a+(b2*a2-2*b*a2*d+d2*a2+4*c*b*e*a)(1/2), b/a 1/2/c/b*(-b*a+d*a-(b2*a2-2*b*a2*d+d2*a2+4*c*b*e*a)(1/2), b/aplot(z(1,:),z(2,:),equilibrium(1,1),equilibrium(1,2),r*)title(系统演化曲线: 相空间)xlabel(Numbers of Patient)ylabel(Numbers of Easy to be Infected)附应用程序Runge_Kutta.mfunction varargout = Runge_Kutta(t,y,flag,a,b,c,d,e)switch flagcase % Return dy/dt = f(t,y). varargout1 = f(t,y,a,b,c,d,e);case init % Return default tspan,y0,options. varargout1:3 = init;case jacobian % Return Jacobian matrix df/dy. varargout1 = jacobian(t,y,a,b,c,d,e);otherwise error(Unknown flag flag .);end% -function dydt = f(t,y,a,b,c,d,e)dydt =a*y(1,:).*y(2,:)-b*y(1,:) -a*y(1,:).*y(2,:)-c*(y(1,:).2).*y(2,: