微分方程数值解--大纲

第十章 偏微分方程数值解法偏微分方程问题，其求解十分困难。

除少数特殊情况外，绝大多数情况均难以求出精确解。

因此，近似解法就显得更为重要。

本章仅介绍求解各类典型偏微分方程定解问题的差分方法。

§1 差分方法的基本概念1.1 几类偏微分方程的定解问题椭圆型方程：其最典型、最简单的形式是泊松（Poisson ）方程),(2222y x f yu x u u =∂∂+∂∂=∆ 特别地，当0),(≡y x f 时，即为拉普拉斯（Laplace ）方程，又称为调和方程2222=∂∂+∂∂=∆yux u u Poisson 方程的第一边值问题为⎪⎩⎪⎨⎧Ω∂=Γ=Ω∈=∂∂+∂∂Γ∈),(),(),(),(),(2222y x y x u y x y x f y ux u y x ϕ 其中Ω为以Γ为边界的有界区域，Γ为分段光滑曲线，ΓΩ称为定解区域，),(y x f ，),(y x ϕ分别为Ω，Γ上的已知连续函数。

第二类和第三类边界条件可统一表示为),(),(y x u u y x ϕα=⎪⎪⎭⎫ ⎝⎛+∂∂Γ∈n 其中n 为边界Γ的外法线方向。

当0=α时为第二类边界条件， 0≠α时为第三类边界条件。

抛物型方程：其最简单的形式为一维热传导方程220(0)u ua a t x∂∂-=>∂∂ 方程可以有两种不同类型的定解问题：初值问题⎪⎩⎪⎨⎧+∞<<∞-=+∞<<-∞>=∂∂-∂∂x x x u x t x u a tu )()0,(,0022ϕ初边值问题221200,0(,0)()0(0,)(),(,)()0u ua t T x l t x u x x x lu t g t u l t g t t Tϕ⎧∂∂-=<<<<⎪∂∂⎪⎪=≤≤⎨⎪==≤≤⎪⎪⎩其中)(x ϕ，)(1t g ，)(2t g 为已知函数，且满足连接条件)0()(),0()0(21g l g ==ϕϕ边界条件)(),(),(),0(21t g t l u t g t u ==称为第一类边界条件。

随机微分方程的数值解引言随机微分方程（Stochastic Differential Equation，简称SDE）是描述包含随机变量的微分方程，它在金融、物理学、生物学等领域具有广泛的应用。

与确定性微分方程相比，SDE中的随机项引入了不确定性和随机性，使得问题更具挑战性和现实性。

本文将介绍随机微分方程的基本概念、求解方法和数值解的计算。

一、随机微分方程概述1.1 确定性微分方程与随机微分方程的区别•确定性微分方程：一般形式为 dy(t) = f(y(t), t)dt，其中f是已知的函数，表示因变量y的增量与自变量t的关系。

•随机微分方程：一般形式为 dy(t) = f(y(t), t)dt + g(y(t), t)dW(t)，其中dW(t)是一个随机项，通常表示为Wiener过程或布朗运动。

1.2 随机微分方程的数学表达一般形式的随机微分方程可以表示为： dy(t) = f(y(t), t)dt + g(y(t),t)dW(t)，其中： - y(t)是待求解的随机过程； - f(y(t), t)表示因变量y的增量与自变量t之间的确定性关系； - g(y(t), t)表示因变量y的增量与自变量t 之间的随机关系； - dW(t)是一个随机项，通常表示为Wiener过程或布朗运动。

二、随机微分方程的求解方法2.1 解析解方法对于简单形式的随机微分方程，可以通过解析的方法求得解析解。

然而，大多数情况下，由于随机视频和随机关系的存在，解析解并不存在或难以求得。

2.2 数值解方法数值解是求解随机微分方程的主要方法之一，它通过将时间间隔分割为若干小段，采用数值方法近似求解微分方程。

常用的数值解方法有： 1. 欧拉方法（Euler Method）：将时间间隔分割为若干小段，在每个小段内使用线性逼近的方式求解微分方程。

2. 随机插值方法（Stochastic Interpolation Method）：利用数值差分逼近计算随机项的变化，并采用插值方法求解微分方程。

第十章 偏微分方程数值解法一、 典型的偏微分方程介绍 1．椭圆型方程 科学技术中经常遇到一些重要的、典型的偏微分方程。

在研究有热源稳定状态下的热传导，有固定外力作用下薄膜的平衡问题时，都会遇到Poisson 方程D y x y x f yux u ∈=∂∂+∂∂),(),(2222（10.1）其中D 表示平面区域。

特别在没有热源或没有外力时，就得到Laplace 方程02222=∂∂+∂∂y ux u （10.2）此外，当研究不可压缩理想流体无旋流动的速度势以及静电场的电位等，也会遇到（10.1）或（10.2）类型的方程。

2．抛物型方程 在研究热传导过程、气体扩散现象、电磁场的传播等问题中以及在统计物理、概率论和重子力学中，经常遇到抛物型方程。

这类方程中最简单、最典型的是热传导方程。

L x t xu a t u <<>=∂∂-∂∂0,0,022（10.3）其中a 是常数。

它表示长度为L 的细杆内，物体温度分布的规律。

3．双曲型方程 在研究波的传播、物体的振动时，常遇到双曲型方程。

这类方程中最简单、最典型的是波动方程L x t xu a t u <<>=∂∂-∂∂0,0,022222（10.4）它表示长度为L 的弦振动的规律。

二、定解问题偏微分方程（10.1）~（10.4）是描述物理过程的普遍规律的。

要使它们刻划某一特定的物理过程，必须给出附加条件。

把决定方程唯一解所必须给定的初始条件和边界条件叫做定解条件。

定解条件由实际问题提出。

对方程（10.3）来说，初始条件的提法应为)()0,(x f x u =，其中f (x )为已知函数，它表示物体在初始状态下温度分布是已知的。

边界条件的提法应为物体在端点的温度分布为已知，即⎩⎨⎧≥==0)(),()(),0(t t t L u t t u ψϕ （10.5）其中ϕ（t ）和ψ（t ）为已知函数。

对（10.4）来说，边界条件的提法和（10.5）形式一样，它表示弦在两端振动规律为已知。

中国海洋大学本科生课程大纲课程属性：公共基础/通识教育/学科基础/专业知识/工作技能，课程性质：必修、选修一、课程介绍1.课程描述：本课程介绍数值求解偏微分方程的基本方法及相关的理论基础。

本课程针对数学类专业高年级（三年级）本科生开设。

课程基本内容包括：有限差分方法、差分格式的稳定性、收敛性分析；变分原理，Galerkin有限元方法等。

通过对模型问题的基本数值方法进行分析，阐明构造数值方法的基本思想和技巧。

通过本课程学习，使学生了解并掌握数值求解偏微分方程的基本思想、基本概念和基本理论（数值格式的相容性、稳定性、收敛性及误差估计等），能够运用算法语言对所学数值方法编制程序在计算机上运行实施并对数值结果进行分析。

培养学生理论联系实际，解决实际问题的能力和兴趣。

2.设计思路：偏微分方程是应用数学的核心内容，在其他科学、技术领域具有广泛深入的应用。

掌握偏微分方程的基础理论及求解方法是数学类专业本科生培养的基本要求。

本课程是在数学物理方程课程基础上开设的延展应用型课程，是一门数值分析理论与实践应用高度融合的专业课。

课程引导学生通过数值方法探讨和理解应用数学工具解决实际- 6 -问题的途径及理论分析框架。

学习本课程需要学生掌握了“数学分析”、“数学物理方程”、“数值分析”及“泛函分析”的核心基本内容。

课程内容安排分为有限差分方法和有限元方法两个单元模块，这是目前应用最广泛、理论发展最完善的两类数值方法，两者既有关联又有本质区别，能够体现偏微分方程数值解法的基本特征。

首先介绍有限差分方法。

有限差分方法是近似求解偏微分方程的应用最广泛的数值方法，以对连续的“导数（微分）”进行离散的“差分”近似为基本出发点，利用Fourier 分析及数值分析的基本理论，讨论椭圆、抛物、双曲等三类典型偏微分方程近似求解方法及近似方法的数学理论分析。

有限元方法是20世纪中期发展起来的基于变分原理的数值方法，具有更直接的物理背景含义，因而受到力学、工程等应用领域广泛的关注和应用。

微分方程的解析与数值解法微分方程既是数学分析的重要分支，也是许多学科领域的基础。

在实际问题的求解中，我们常常需要寻找微分方程的解析解或者数值解。

本文将围绕微分方程的解析和数值解法展开讨论。

一、微分方程的解析解解析解指的是通过代数计算得到的方程的解。

对于某些简单的微分方程，我们可以通过分离变量、变量代换等方法得到解析解。

下面以一阶线性常微分方程为例，讨论解的求解过程。

考虑一阶线性常微分方程形式如下：$$\frac{dy}{dx} + P(x)y = Q(x)$$其中，$P(x)$和$Q(x)$为已知函数。

我们可以通过以下步骤求解该微分方程：1. 将方程改写为标准形式：$\frac{dy}{dx} + P(x)y - Q(x) = 0$2. 求解齐次线性微分方程：$\frac{dy}{dx} + P(x)y = 0$。

记其解为$y_h$，即$y_h = Ce^{-\int P(x)dx}$，其中$C$为常数。

3. 利用常数变易法，假设原方程的解为$y = u(x)y_h$，其中$u(x)$为待定函数。

4. 将$y = u(x)y_h$代入原方程，得到关于$u(x)$的方程。

5. 求解$u(x)$的方程，得到$u(x)$的表达式。

6. 将$u(x)$代入$y = u(x)y_h$，得到原方程的解析解。

上述过程就是一阶线性常微分方程求解的一般步骤。

对于其他类型的微分方程，也有相应的解析解求解方法。

但并非所有微分方程都存在解析解。

二、微分方程的数值解法对于一些复杂的微分方程，无法找到解析解，此时我们需要借助数值方法求解。

常见的数值解法包括欧拉法、改进的欧拉法、四阶龙格-库塔法等。

1. 欧拉法欧拉法是一种较为简单的数值解法，其基本思想是通过离散化微分方程，将微分方程转化为差分方程。

具体步骤如下：将求解区间$[a, b]$等分成$n$个小段，步长为$h = \frac{b-a}{n}$。

利用微分方程的导数定义，将微分方程转化为差分方程，即$y_{i+1} = y_i + h \cdot f(x_i, y_i)$，其中$f(x, y)$为微分方程右端的函数。

微分方程的常用数值解法摘要：微分方程是数学中的一种重要的方程类型，它能描述自然现象和工程问题中的许多变化规律。

但是大多数微分方程解法是无法用解析的方式求解的，因此需要借助数值解法来近似求解。

本文将介绍微分方程的常用数值解法。

关键词：欧拉方法；龙格-库塔方法；微分方程；常用数值解法一、微分方程数值解方法微分方程数值解法是数学中的重要部分。

欧拉方法、龙格-库塔方法和二阶龙格-库塔方法是常用的微分方程数值解法，下面就分别介绍这三种方法。

（一）欧拉方法欧拉方法是解初值问题的一种简单方法，它是欧拉用的第一种数值方法，也叫向前欧拉法。

欧拉方法是利用微分方程的定义式y’=f(x, y)，将它带入微分方程初值问题y（x_0）=y_0中，以y_0为初始解，在每一步上通过沿着切线的方法进行估计并推进新的解y_{i+1}：y_i+1=y_i+hf(x_i,y_i)其中，x_i和y_i是我们知道的初始条件，h是求解过程中的步长，f是微分方程右端项。

它是一种时间迭代的算法，易于实现，但存在着精度不高的缺点。

（二）龙格-库塔方法龙格-库塔方法是一种经典迭代方法，也是近代微分方程数值解法发展的里程碑之一。

龙格-库塔方法的主要思想是利用规定的阶码及阶向量，通过递推求解微分方程数值解的近似值。

龙格-库塔方法的方式不同，其步骤如下：第一步：根据微分方程，计算出在x_i和y_i的值。

第二步：在x_i处对斜率进行估计，并利用这个斜率来求解下一步所需的y_i+1值。

第三步：使用x_i和y_i+1的值来重新估计斜率。

第四步：使用这个新的斜率来更新y_i+1的值。

（三）二阶龙格-库塔方法二阶龙格-库塔方法是龙格-库塔方法的一种变体，它根据龙格-库塔方法的思想，使用更好的步长来提高数值解的精度。

二阶龙格-库塔方法的基本思路是，在第一次迭代时使用一个阶段小一半的y_i+1，然后使用这个估算值来计算接下来的斜率。

通过这种方法，可以提高解的精度。

二阶龙格-库塔方法的步骤如下：第一步：计算出初始阶段的y_i+1值。

常微分方程课程教学大纲（Ordinary Differential Equation）课程性质：学科基础课适用专业：信息与计算科学先修课程：数学分析、高等代数、普通物理后续课程：微分方程数值解总学分：3教学目的与要求：微分方程是数学理论联系实际的重要渠道之一，也是其它数学分支的一个综合应用场所，我们所研究的方程多数是由其它学科（如物理、气象、生态学、经济学）推导而来，通过本课程的学习使学生了解到微分方程和其它数学分支的联系及其在其它自然科学学科中的应用，使学生进一步了解到数学的重要性和广泛的应用背景。

通过对微分方程发展史的回顾，让学生从一个侧面了解人类对自然界的认识过程和科学研究的探索过程，逐步培养学生的活学活用能力和创造发展的能力。

通过本课程的学习，使学生熟练掌握各类方程的判别与求解，掌握基本理论的基本思想和证明方法。

并简要介绍一些其它学科需要我们解决而目前我们尚不能解决的问题，为其它后续课程留下引子，并通过一些例子让学生知道目前这个学科的最新研究动态。

教学内容与学时安排第一章绪论（4学时）1、微分方程的产生；2、基本概念。

第二章一阶微分方程的初等解法（12学时）1、变量分离方程与变量变换；2、线性方程与常数变易法；3、恰当方程与积分因子；4、一阶隐方程与参数表示。

本章重点：各种初等解法。

难点：一阶隐方程与参数表示。

第三章一阶微分方程解的存在定理（6学时）1、解的存在唯一性定理与逐步逼近法；2、解的延拓。

本章重点：解的存在唯一性定理与逐步逼近法。

难点：解的延拓。

第四章高阶微分方程（10学时）1、线性微分方程的一般理论；2、常系数线性方程的解法；3、高阶方程的降阶和幂级数解法。

本单重点：常系数线性方程的解法。

难点：高阶方程幂级数解法。

第五章线性微分方程组（10学时）1、存在唯一性定理；2、线性微分方程组的一般理论；3、常系数线性微分方程组。

本章重点：线性微分方程组的一般理论。

第六章非线性微分方程和稳定性（10学时）1、引言；2、相平面；3、按线性近似决定微分方程组的稳定性；4、李雅普诺夫第二方法；5、周期解和极限环。

微分方程数值解法微分方程数值解法是一种将微分方程的解转化为数值计算的方法。

常用的微分方程数值解法包括欧拉法、隐式欧拉法、龙格-库塔法等。

1. 欧拉法：欧拉法是最简单的一种数值解法，它基于微分方程的定义，在给定的初始条件下，通过不断迭代计算微分方程在给定区间上的近似解。

欧拉法的迭代公式为：y_{n+1}=y_n+h\\cdot f(t_n,y_n)，其中y_n表示第n步的近似解，t_n表示第n步的时间，h表示步长，f(t_n,y_n)表示微分方程的右侧函数。

2. 隐式欧拉法：隐式欧拉法是欧拉法的改进，它在计算近似解时使用了未知公式的近似值，从而提高了精度。

隐式欧拉法的迭代公式为：y_{n+1}=y_n+h\\cdotf(t_{n+1},y_{n+1})，其中y_{n+1}表示第n+1步的近似解，t_{n+1}表示第n+1步的时间，h表示步长，f(t_{n+1},y_{n+1})表示微分方程的右侧函数。

3. 龙格-库塔法：龙格-库塔法是一种常用的高阶数值解法，它通过计算微分方程的斜率来提高精度。

最常见的是四阶龙格-库塔法，它的迭代公式为：y_{n+1}=y_n+\\frac{1}{6}(k_1+2k_2+2k_3+k_4)，其中k_1=h\\cdot f(t_n,y_n)，k_2=h\\cdotf(t_n+\\frac{h}{2},y_n+\\frac{1}{2}k_1)，k_3=h\\cdotf(t_n+\\frac{h}{2},y_n+\\frac{1}{2}k_2)，k_4=h\\cdotf(t_n+h,y_n+k_3)。

这些方法的选择取决于问题的性质和精度要求。

其中，欧拉法是最简单的方法，但精度较低，龙格-库塔法精度较高，但计算量较大。

在实际应用中需要根据问题的具体情况选择合适的数值解法。

第一章绪论1.1 引言常微分方程是现代数学的一个重要分支，是人们解决各种实际问题的有效工具。

微分方程的理论和方法从17世纪末开始发展起来，很快成了研究自然现象的强有力工具，在17到18世纪，在力学、天文、科学技术、物理中，就已借助微分方程取得了巨大的成就。

1864年Leverrer根据这个方程预见了海王星的存在，并确定出海王星在天空中的位置。

现在，常微分方程在许多方面获得了日新月异的应用。

这些应用也为常微分方程的进一步发展提供了新的问题，促使人们对微分方程进行更深入的研究，以便适应科学技术飞速发展的需要。

研究常微分方程常用数值解是数学工作者的一项基本的且重要的工作。

在国内外众多数学家的不懈努力，使此学科基本上形成了一套完美的体系。

微分方程的首要问题是如何求一个给定方程的通解或特解。

到目前为止，人们已经对许多微分方程得出了求解的一般方法。

由于在生产实际和科学研究中所遇到的微分方程问题比较复杂，使这些问题的解即使能求出解析表达式，也往往因计算量太大而难于求出，而对于一些典型的微分方程则可以运用基本方法求出其解析解，并可以根据初值问题的条件把其中的任意常数确定下来。

由于求通解存在许多困难，人们就开始研究带某种定解条件的特解。

首先是Cauchy对微分方程初始解的存在惟一性进行了研究。

目前解的存在惟一性、延拓性、大范围的存在性以及解对初始解和参数的延续性和可微性等理论问题都已发展成熟。

与此同时，人们开始采取各种近似方法来求微分方程的特解，例如求微分方程数值解的Euler折线法、Runge-Kutta法等，可以求得若干个点上微分方程的近似解。

最后，由于当代高科技的发展为数学的广泛应用和深入研究提供了更好的手段。

用计算机结合Matlab软件求方程的精确解、近似解，对解的性态进行图示和定性、稳定性研究都十分方便有效。

本章先介绍常微分的一般概念、导出微分方程的一些典型例子及求解微分方程的思路分析。

从而得到常微分方程的常用数值解法。

微分方程数值解微分方程数值解及其应用绪论自然界中的许多事物的运动和变化规律都可以用微分方程来描述，因此对工程和科学技术中的实际问题的研究中, 常常需要求解微分方程．但往往只有少数较简单和典型的微分方程可求出其解析解，在大多数情况下,只能用近似法求解，数值解法是一类重要的近似方法．本文主要讨论一阶常微分方程的初值问题的数值解法，探讨这些算法在处理来自生活实际问题中的应用，并结合MATLAB 软件，动手编程予以解决．１ 微分方程的初值问题[1] 1.1 预备知识在对生活实际问题的研究中，通常需要考虑一阶微分方程的初值问题00(,)()dyf x y dx y x y ⎧=⎪⎨⎪=⎩（1）这里(),f x y 是矩形区域R ：00,x x a y y b -≤-≤上的连续函数．对初值问题（1）需要考虑以下问题：方程是否一定有解呢?若有解，有多少个解呢?下面给出相关的概念与定理．定义 1 Lipschitz 条件[1][2]：矩形区域R ：00,x x a y y b -≤-≤上的连续函数(),f x y 若满足：存在常数0L >，使得不等式()()1212,,f x y f x y L y y -≤-对所有()()12,,,x y x y R ∈都成立，则称(),f x y 在R 上关于y 满足Lipschitz 条件.定理 1 解的存在唯一性定理[1][3]：设f 在区域()}{,,D x y a x b y R =≤≤∈上连续，关于y 满足Lipschitz 条件，则对任意的[]00,,∈∈x a b y R ，常微分方程初值问题（1）当[],x a b ∈时存在唯一的连续解()y x ．该定理保证若一个函数(),f x y 关于y 满足Lipschitz 条件，它所对应的微分方程的初值问题就有唯一解.在解的存在唯一性得到保证的前提下，自然要考虑方程的求1()()()(,)n n n n n y x y x y x f x y h+-'≈=将初值问题（1）离散化，则问题（1）可化为1(,),n n n n y y h f x y +=+⋅0n x x n h =+⋅，（3）此方法称为欧拉法.欧拉方法的几何意义在数值计算公式中体现了出来.在xy 平面上，一阶微分方程的解()y y x =称作它的积分曲线.积分曲线上一点(),x y 的切线斜率等于函数(),f x y ，按函数(),f x y 在xy 平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致.基于上述几何解释，从初始点000(,)P x y 出发，先依方向场在该点的方向上推进到1x x =上一点1P ，再从1P 依方向场的方向推进到2x x =上一点2P ，循环前进便作出一条折线012P PP ，因此欧拉方法又称为折线法.若初值0y 已知，则由（3）式可逐步算出10002111(,),(,),y y h f x y y y h f x y =+⋅=+⋅为了分析计算公式的精确度，通常可用泰勒展开将()1n y x +在n x 处展开，则有()()()()()()2''11,,.2++'=+=++∈n n n n n n n n h y x y x h y x y x h y x x ξξ在()n n y y x =的前提下，()()()(),,.n n n n n f x y f x y x y x '==可得欧拉法（3）的误差为()()()2211.22n n n n h h y x y y y x ξ++''''-=≈容易看出，欧拉法（3）式具有一阶精度.2）向后欧拉方法：如果对微分方程（1）从n x 到1n x +积分，得 ()()()()11,n nx n n x y x y x f t y t dx ++=+⎰，（4）如果（4）式右端积分用右矩形公式()()11,n n h f x y x ++⋅近似，则得到另一个公式 ()111,n n n n y y hf x y +++=+， （5） 称为后退欧拉法.值得一提的是:后退欧拉法与欧拉公式有着本质的区别，后者是关于1n y +的直接计算公式，它是显式的，而（5）式的右端含有关于1n y +的表达式，它是隐式的.在利用后退欧拉法时，我们通常利用迭代法求解，实质就是逐步显示化.具体迭代过程如下：首先利用欧拉公式(0)1(,)+=+⋅n n n n y y h f x y 给出迭代初值(0)1+n y ，把它代入（5）式的右端，使之转化为显式，直接计算得11(1)(0)1(,)+++=+⋅n n n n y y h f x y .如此反复进行，得 (1)()111(,)++++=+⋅k k n n n n y y h f x y 0,1,k = ，则得到后退欧拉法的迭代公式(0)1(1)()111(,)(,)+++++⎧=+⋅⎨=+⋅⎩nn n n k k n n n n y y h f x y y y h f x y ， 可以看出，后退欧拉法具有一阶精度，且计算比较麻烦. 1.2.2梯形方法为得到比欧拉法精确度高的计算公式，在等式（4）式右端积分中若用梯形求积公式近似，并用n y 代替()n y x ，1n y +代替()1n y x +，则得()()111,,2n n n n n n hy y f x y f x y +++=++⎡⎤⎣⎦， （6）称其为梯形方法.梯形方法与后退欧拉法一样，都是隐式单步法，可用迭代法求解，其迭代公式为()()(0)1(1)()111(,),,2+++++⎧=+⋅⎪⎨⎡⎤=++⎪⎣⎦⎩n n n n k k n n n n n n y y h f x y h y y f x y f x y . （7） 为了分析梯形公式的收敛性，将（6）与（7）式相减，得()()(1)()111111,,2k k n n n n n n h y y f x y f x y +++++++⎡⎤-=-⎣⎦，0,1,2,k= 因为(),f x y 满足Lipschitz 条件，于是有(1)()11112+++++-≤-k k n n n n hL y y y y ，其中L 为(),f x y 关于y 的Lipschitz 常数.如果选取h 充分小，使得12hL<，则当k →∞时有(1)11+++→k n n y y ，这说明迭代过程（7）式是收敛的[4].容易推导得出梯形法（7）式是二阶方法.经分析，梯形方法虽然提高了精度，但是以增加计算量为代价的.从上述的迭代公式可以看出，每迭代一次都要重新计算(),f x y 的值，而且迭代又要进行若干次，计算相当的复杂.为此，有没有比较简便的计算方法呢？下面给出改进的欧拉方法. 1.2.3改进的欧拉方法由前面的讨论可知，梯形法计算相对复杂，现对上面的梯形法进行简化，具体方法是只计算一两次就转入下一步的计算，先用欧拉公式（3）求得一个初步的近似解1n y +，称为预测值，再利用公式（6）把它校正一次，这样建立的预测-校正系统通常称为改进的欧拉公式.具体公式如下()()()11,,,2n n n n n n n n hy y f x y f x y hf x y ++⎡⎤=+++⎣⎦（8）改进的欧拉法与梯形法一样，是二阶方法. 1.2.4 Runge-Kutta 方法由前面讨论可知，从（4）式()()()()11,n nx n n x y x y x f x y x dx ++-=⎰可以看出，若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，它必然要增加求积积点，为此将（4）式的右端用求积公式表示为()()()()11,,n nrx i n i n i x i f x y x dx h c f x h y x h λλ+=≈++∑⎰，（9）一般来说，点数r 越多，精度越高，上式右端相当于增量函数(),,x y h ϕ，为得到便于计算的显式方法，将公式（9）表示为：()1,,,n n n n y y h x y h ϕ+=+（10） 其中()()1111,,,,,2,r n n i i i n n i i n i n ij j j x y h c K K f x y K f x h y h K i rϕλμ=-=⎧⎪=⎪⎪=⎨⎪⎛⎫⎪=++= ⎪⎪⎝⎭⎩∑∑ （11）这里,,i i ij c λμ均为常数. i c 为加权因子，i K 为第i 段斜率，共有r 段.我们把（10）和（11）称为r 级显式Runge-Kutta 法，简称为R-K 方法.下面给出其中最经典最常用的一个公式：()()()11234121324322,6,,,22,,22,.n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎛⎫⎪=++ ⎪⎨⎝⎭⎪⎪⎛⎫=++ ⎪⎪⎝⎭⎪⎪=++⎩（12）Runge-Kutta 方法作为一种重要的单步方法，具有很高的实用价值，它关于初值是稳定的，其解连续地依赖于初值，是一类便于应用的单步法，为了计算1n y +，只用到前面一步的值ny 即可，因此每步的步长可以独立取定.常用的Runge-Kutta 方法精度较高，为了达到预定的精度，与欧拉方法与梯形法相比，步长h 可取得大些，求解区间上的总步数可以少些.但Runge-Kutta 方法也有些缺点，比如四阶Runge-Kutta 方法每算一步需要四次计算(),f x y 的值，计算量较大（对于复杂的(),f x y 而言）．2 数值方法的应用实例[5-9]例1 对于初值问题()1001y yy '=-⎧⎪⎨=⎪⎩，分别用欧拉法、改进的欧拉法，梯形法求()1y 的近似值.解：易得该方程的解析解()10x=，()1 4.5400e-005y=，为比较，将按不同y x e-数值计算方法所得结果列表如下：表 1 三种不同方法的数值结果0.20.40.60.8100.0050.010.0150.02图 1 三种不同方法数值解与精确解的误差曲线从表1中可以看出：当0.2h =时，三种方法均不稳定，计算结果严重偏离精确值；0.1h =时，改进后的欧拉和梯形法均稳定，但欧拉法效果很差；当0.01h ≤时，三种方法均稳定，但精确度有区别．可以看出，h 越小，计算结果越好，要想计算结果充分接近于解析解还须取较小的h 值．图1反映的步长0.01h =时，三种数值方法的所得数值解与解析解在[0,1]区间的误差曲线，由图可知，在步长相同的情况下，梯形法的精确度略高于改进的欧拉法；改进的欧拉法和梯形法精确度都明显高于欧拉法．例2 用欧拉法、改进的欧拉法和Runge-Kutta 法求解初值问题[]()2,0,101dyx y x dxy y ⎧=-∈⎪⎨⎪=⎩并比较三种方法的结果．解：方程为1n =-的伯努利方程，可求得解析解为2y e =现用MATLAB 软件编程，用题目要求的方法求解，可得如下图示结果：图2 （a）步长为0.2时R-K 法和解析解比较图2 （b）步长为0.2时改进的Euler法和解析解比较图2 （c）步长为0.2时欧拉法和解析解比较上图2(a)，(b)，(c)描述的是步长为0.2时，用欧拉法、改进的欧拉法，Runge-Kutta法求解方程所得的数值解与解析解之间的对比图．由图可知，Runge-Kutta法所得数值解曲线和解析解曲线吻合的很好，改进的欧拉法和欧拉法随着计算的进行，数值解和解析解之间误差逐步增大，但改进的欧拉法效果要好于欧拉法．图3 (a) 步长为0.1时Euler法和解析解比较图3 (b) 步长为0.1时改进的Euler法和解析解比较图3 (c) 步长为0.1时Runge-Kutta法和解析解比较上图3 (a)，(b)，(c)描述的是步长为0.1时，用欧拉法、改进的欧拉法，Runge-Kutta 法求解方程所得的数值解与解析解之间的对比图．由图可知，改进的欧拉法和Runge-Kutta法所得数值解曲线和解析解曲线吻合的很好，而欧拉法随着计算的进行数值解和解析解之间误差逐步增大．相应的程序如下：主程序x=0:0.2:1;jxj=exp(2*x).*(1./exp(4*x) + (2*x)./exp(4*x)).^(1/2);y=Euler(@ff,0,1,0.2,1);gy=geuler(@ff,0,1,0.2,1);Ry=RK(@ff,0,1,0.2,1);figure(1);plot(x,jxj,x,Ry,'*');figure(2);plot(x,jxj,x,gy,'*');figure(3);plot(x,jxj,x,y,'*')欧拉法程序function y=Euler(f,a,b,h,y0)n=(b-a)/h;x=a:h:b;y=zeros(n+1,1);y(1)=y0;for i=1:ny(i+1)=y(i)+h*feval(f,x(i),y(i));end改进的欧拉法程序function gy=geuler(f,a,b,h,y0)n=(b-a)/h;x=a:h:b;y=zeros(n+1,1);y(1)=y0;for i=1:nyp=y(i)+h*feval(f,x(i),y(i));yc=y(i)+h*feval(f,x(i+1),yp);y(i+1)=(yp+yc)/2;endgy=y;Runge-Kutta法程序function Ry=RK(f,a,b,h,y0)n=(b-a)/h;x=a:h:b;y=zeros(n+1,1);y(1)=y0;for i=1:nk1=feval(f,x(i),y(i));k2=feval(f,x(i)+h/2,y(i)+h*k1/2);k3=feval(f,x(i)+h/2,y(i)+h*k2/2);k4=feval(f,x(i+1),y(i)+h*k3);y(i+1)=y(i)+h*(k1+2*k2+2*k3+k4)/6;endRy=y;3 微分方程数值解法在实际生活中的应用3.1应用实例:耐用消费新产品的销售规律模型一种新产品进入市场以后，常常会经历销售量首先慢慢增加然后逐渐慢慢下降的一个过程，由此给出的时间—销售坐标系下的曲线称为产品的生命曲线,它的形状是钟形的.不过对于较耐用的消费品，情况会有所不一样，它的生命曲线会在开始有个小小的高峰，之后是段比较平坦的曲线，先下降，之后再上升，然后达到高峰，因此它是双峰型的曲线.如何理解这种和传统产品的生命曲线理论相冲突的现象？澳大利亚学者斯蒂芬斯与莫赛观察到的购买耐用消费产品的人大概可分为两类：其一是非常善于接受新的事物，称其为“创新型”消费者，他们会经常从产品广告，制造商给出产品的说明书与商店样品中了解了产品功能与性能之后，再决定其否购买；其二是消费者相对保守，称其为“模仿型”消费者，他们往往会根据大部分已经购买产品的消费者实际使用的经验而提供的信息来决定是否购买产品，下面的例子经过分析，建立相应的数学模型，对这种现象给出了科学解释.3.1.1模型的建立将顾客获得信息大致可分成两类，其一称之为“搜集型”，源自于产告说明、广告，“创新型”顾客获得这类消息后就可作出其是否购买；第二类称之为“体验型”，即消费者使用之后会获得真实体验，常常以口头的形式散播，“模仿型”类顾客会在获得这种信息之后才能决定是否购买.设K 是潜在用户的总数，1K 与2K 分别为 “创新型”与“模仿型”的人数.再设()N t 是时刻t 已经购买的商品顾客的人数，而()1N t 与()2N t 分别为“创新型”与“模仿型”的顾客的人数,再设()1A t 是时刻t 中已获得“搜集型”信息的人数，由于该部分的信息可直接由外部获得，或者可从已获得该类信息的人群中获得，因此类似巴斯模型，从而建立如下方程：()()()()()()11112112,00,,0dA t K A t A t A dtαααα=-+=> ， 获得“搜集型”信息的“创新型”消费者决定其是否购买的行为，满足如下方程：()()()()()()11111,00,,0dN t K N t N t N dtαβαβ=-+=> 对于“模仿型”的顾客，可从已经购买该产品“创新型”或者“模仿型”的顾客中获取信息，于是有()()()()()()22212,0dN t K N t N t N t dtγγ=-+> 在这里，忽略消费者购买产品后需一段短暂的使用后才会散播体验信息的滞后作用． 综上，斯蒂芬斯—莫赛模型是一常微分方程组的初值问题：()()()()()()()()()()(()11112221212,00,(0)0dN t K N t N t dt dN t K N t N t N t dt N N αβγ⎧=-+⎪⎪⎪=-+⎨⎪==⎪⎪⎩而()()()12N t N t N t =+为时刻t 购买该商品的总人数. 3.1.2 模型的求解对于斯蒂芬斯—莫赛模型中()2N t 的解析解则不能求出，于是可以用四阶Runge Kutta -公式求得，且在它的精度要求达到很高情形下求出()2N t .用MATLAB 软件求解，相应的程序及结果如下function RK=RKFC(fc,a,b,h,y0)n=(b-a)/h;x=a:h:b;m=length(y0);y=zeros(n+1,m);y(1,:)=y0; for i=1:nk1=feval(fc,x(i),y(i,:));k2=feval(fc,x(i)+h/2,y(i,:)+h*k1/2); k3=feval(fc,x(i)+h/2,y(i,:)+h*k2/2);k4=feval(fc,x(i+1),y(i,:)+h*k3); y(i+1,:)=y(i,:)+h*(k1+2*k2+2*k3+k4)/6; end RK=y;function f=FC(x,y)k1=50; k2=70;al=0.0013;be=0.0013;ga=0.0015;f=[(k1-y(1))*(al+be*y(1)),ga*(k2-y(2))*(y(1)+y(2))];x=0:0.3:24;RK=RKFC(@FC,0,24,0.3,[0,0])figure(1)；plot(x,RK (:,1),'+',x,RK(:,2),'*')；legend('N1(t)','N2(t)',2) figure(2)；plot(x,RK(:,1)+RK(:,2),'+')；legend('N(t)',2)图4()()12,N t N t 与时间关系图图 5 ()N t 与[0,25]时间段关系图由此例可以看出，微分方程数值解在实际生活有着广泛的应用，而数值解法中的Runge-Kutta方法是一种很重要且应用很广泛的算法.结语微分方程数值解是求解微分方程的一种很重要且应用范围很广的方法，而常用的数值解法如欧拉法、改进的欧拉法、梯形法和Runge-Kutta方法在一定程度上都有自己的优缺点，理解和掌握各种方法的应用范围，用MATLAB编制各种方法求解实际问题的通用程序，对用微分方程数值解理论解决现实生活中的实际问题有很重要的意义．参考文献[1] 李庆扬,王能超,易大义.数值分析[M].北京:清华大学出版社,清华大学出版社,2001.[2] 冯康.数值计算方法[M].杭州:浙江大学出版社,2003.[3] 封建湖.计算方法典型题分析解集[M].西安:西北工业大学出版社,2000.[4] 胡建伟,汤怀民.微分方程数值解法(第二版)[M].北京:科学出版社,2007.[5] 王能超.计算方法简明教程[M].北京:高等教育出版社,2004.[6] Nash S G.A history of scientific computing.[M] New York:ACM Press,1990.[7] 于丽妮. ODE问题解析解及数值解的matlab实现[J]. 电脑知识与技术.2012,8(14):3303-3305.[8] 霍晓成. 常微分方程数值解法的研究[J]. 临沂师范学院学报. 2011,(6):19-23.[6] 王国立, 陈瑛.非线性微分方程迭代算法及其在物理学中的应用[J]. 长春师范学院学报( 自然科学版).2006, 25(2):10-12.。