[常微分方程模型及其数值解2008.ppt

常微分方程数值解 杨瑞琰 0、导　　言 在许多实际问题中，例如物理中的速率问题，人口的增长问题，放射性衰变问题，经济学中的边际问题等，常常涉及到两个变量之间的变化规律。微分方程是研究上述问题的一种机理分析方法，它在科技、工程、生态、环境、人口以及经济管理等领域中有着十分广泛的应用。 在应用微分方程解决实际问题时，必须经过两个阶段。一是微分方程的建立，建立一个微分方程的实质就是构建函数、自变量以及函数对自变量的导数之间的一种平衡关系。而正确地构建这种平衡关系，需要对实际问题的深入浅出的刻画，根据物理的和非物理的原理、定律或定理，作出合理的假设和简化并将它升华成数学问题。 另一个是方程的求解和结果分析。对一些常系数的或特殊函数形式的微分方程，往往能得到解析解，这对实际问题的分析和应用都是有利的,但是大多数变系数的、非线性函数形式的微分方程都是求不出解析解的,此时就需要应用求解微分方程的另一个重要方法──数值解法。 本章简要介绍有关微分方程模型的概念，微分方程的数值解法和图解法，主要介绍若干建模实例，通过它们展示微分方程模型的建模步骤及解决实际问题的全过程。 欧拉方法和龙格—库塔方法 function ex03() %用三种不同方法解微分方程 clear, clf x0=0; xt=1; Num=100; h=(xt-x0)/(Num-1); x =x0+[0:Num]*h; a = 1; yt = sqrt(1+2*x); %真值解 y0 = 1; %设定函数初值 PointNum = 4; %设定取点数 [x1,ye] = MyEuler(myfun012,x0,xt,y0,PointNum); [x2,yh] = MyEulerPro(myfun012,x0,xt,y0,PointNum); [x3,yr] = MyRunge_Kutta(myfun012,x0,xt,y0,PointNum); plot(x,yt,k, x1,ye,b:, x2,yh,g:, x3,yr,r:) legend(真值,简单欧拉法解,改进欧拉法解,Runge-Kutta法解) hold on plot(x1,ye,bo,x2,yh,b+, x3,yr,r*) R—K方法的MATLAB实现 对于微分方程（组）的初值问题 龙格—库塔方法可用如下MATLAB命令实现其计算： [t，x]=ode23(@f，ts，x0，options) [t，x]=ode45(@f，ts，x0，options) 其中ode23用的是3级2阶龙格—库塔公式，ode45用的是以Runge-Kutta-Fehberg命名的 5级4阶公式。 命令的输入f是待解方程写成的函数m文件： function dx=f(t,x) dx=[f1;f2;?;fn]; 若ts=[t0,t1,t2, …,tf]，则输出在指定时刻t0,t1,t2, …,tf的函数值；等分点时用ts=t0:k:tf，输出在[t0,tf]内等分点处的函数值。x0为函数初值（n维向量）。options可用于设定误差限（options缺省时设定相对误差10-3，绝对误差10-6），命令为： options=odeset(‘reltol’,rt,’abstol’,at) 其中rt，at分别为设定的相对误差和绝对误差。 命令的输出t为指定的ts，x为相应的函数值（n维向量）。 注意，计算步长是根据误差限自动调整的，并不是输入中指定的ts的分点。 由于 yn+1=yn+h?1?n+h?2(?n+?h?xn+?h?n ?yn)+O(h3) O(h3),称此公式为p阶Runge-kutta方法,简称p阶R-K方法. 对于p=2的情形, 应有 =yn+h(?1+?2)?n+h2?2(??xn+??n ?yn)+O(h3) 所以,只要令 ?1+?2=1, ??2=1/2, ??2=1/2 (8) Evaluation only. Created with Aspose.Slides for .NET 3.5 Client Profile 5.2.0.0. Copyright 2004-2011 Aspose Pty Ltd. 一般地, 参数由(8)确定的一族差分公式(7)统称为二阶R-K方法. 称之为中点公式,或可写为 若取?=1,则得?1=?2=1/2,?=1,此时公式(7)就是改进的Euler公式; 若取?1=0,则得?2=1,?=?=1/2,公式(7)为 高阶R-K公式可类似推导. 下面列出常用的三阶、四阶R-K公式. Evaluation only. Crea