偏微分方程的数值方法Read.docx

偏微分方程的数值方法Read.docx

《偏微分方程的数值方法Read.docx》由会员分享，可在线阅读，更多相关《偏微分方程的数值方法Read.docx（13页珍藏版）》请在冰豆网上搜索。

偏微分方程的数值方法Read

第4章偏微分方程的数值方法

在自然科学和工程技术中的许多物理问题，都可以用微分方程来描述。

虽然大学的高等数学课程里研究了一些特殊常微分方程的解析求解，但是在实际应用中大量的微分方程是无法确定其解析解的。

本章和下一章将介绍常微分方程的数值方法、偏微分方程数值求解的思想。

§3.1微分方程数值方法的有关概念

首先介绍微分方程的定义与分类：

含有自变量、未知函数及其导数（微分或偏导数）的方程称为微分方程；

如果未知函数只含有一个变量，则称为常微分方程；如果未知函数含有若干个变量，则称为偏微分方程。

微分方程中未知函数的导数或偏导数的最高阶次称为微分方程的阶。

例如：

微分方程

（3.1.1）

是一阶常微分方程，而

（3.1.2）

是二阶偏微分方程。

所有使微分方程成为等式的函数，都是微分方程的解；在n阶微分方程中，将微分方程的具有n个任意常数的解称为该微分方程的通解。

为确定微分方程通解中的任意常数而列出的条件称为定解条件；定解条件可以分为初始条件和边界条件两类。

由微分方程和定解条件一起构成的问题称为微分方程定解问题。

根据定解条件的不同，常微分方程分为初值问题和边值问题；若定解条件是描述函数在

一点（或初始点）处状态的，则称为初值问题，一阶常微分方程初值问题的一般形式为：

（3.1.3）

若定解条件描述了函数在至少两点（或边界）处状态的称为边值问题，例如：

（3.1.4）

微分方程的解有解析解与数值解两种。

由于实际问题中大多数微分方程无法求出解析解，或难以其用初等函数表示解析解；同

时在应用过程中，一般只需要得到若干特定点处的函数值。

因此微分方程数值解法是科学与

工程计算中的一个重要内容。

由于微分方程数值解是一组离散点处的未知函数值，所以求数值解的基本步骤为：

首先将整个定义域分成若干小块，以便对每小块上的点或片求出近似值，这样按一定规

律对定义域分割的过程称为区域剖分。

其次根据微分方程的形式，构造关于上述离散点或片的函数值递推公式或方程，该步骤

称为微分方程的离散。

这样未知量不再是一个连续函数，而是由若干个未知函数值所构成。

微分方程离散后得到的递推关系式，需要给定若干个初值才能启动。

如果递推式是一个

线性方程组，一般它所含的方程个数要少于未知量的个数，必须补充若干个方程后才可求

解。

这些方程可以通过将微分方程的初始条件或边界条件离散后获得，这一过程称为初始或

边界条件的离散。

经过上面的三个离散化过程，原来的微分方程定解问题就变为离散系统的求解问题。

在

求解之前需要讨论离散系统解的存在唯一性问题；离散系统与微分方程问题之间的差异，即

解的收敛性问题；还需要研究解的收敛速度和计算的稳定性等问题。

最后进行实际计算，通过求解离散系统问题，得到微分方程定解问题的数值解。

应用数学方法解决实际问题可以分为两个阶段，一是对实际问题进行分析，假设，并建

立数学模型；二是根据数学模型的特点，选择适当的数值方法，确定模型的解。

数值方法的

误差主要有以下四个来源。

①模型误差将实际问题归结为数学模型时，需要对问题作一定的简化和假设，由此

产生数学模型与实际问题之间的误差；

②观测误差数学模型中的一些系数、初值等常数源于测量仪器或统计资料，由于客

观条件和仪器精度的限制而产生的误差；

③截断误差数学模型离散化时往往舍去一些次要的项，这将导致数学模型解与离散

问题解之间产生的误差；

④舍人误差利用计算机根据结定的数值方法求解离散问题时，由于计算机对所运算

的对象按一定字长进行四舍五人，将导致问题数值解与离散问题解之间的误差。

本章主要关注微分方程离散化过程中产生的截断误差。

§3.2初值问题的数值方法

本节讨论求解常微分方程（组）初值问题（3.2.1）的数值方法

（3.2.1）

先研究最简单最直观的求解初值问题的Euler方法，然后介绍两类更有效的方法：

Runge－Kutta方法和线性多步方法。

§3.2.1Euler法

由于微分方程的数值解只需要计算在N+1个节点

处微分方程解的近似值．所以先对初值问题（3.2.1）的求解区间

进行剖分，以得到计算节点。

一般将求解区间

均匀分成N等份，即得到的节点

满足：

（3.2.2）

Euler方法的具体计算公式，可以由三种不同的方法推导得到。

（1）差商近似方法

将初值问题（3.2.1）在节点

处的导数

用前向差商代替：

（3.2.3）

记

，则微分方程（3.2.1）近似写成：

（3.2.4）

由初始条件

出发，逐步计算得到

。

式（3.2.4）称为显式Euler公式，显式Euler公式是最基本的计算公式。

如果将节点

处的导数

用后向差商代替：

（3.2.5）

则类似可得递推公式：

（3.2.6）

由于它关于

是隐式形式，所以式（3.2.6）称为隐式Euler公式。

显式和隐式Euler公

式在计算

时，只用到前一步的结果

，可称为单步方法。

如果将节点

处的导数

用中心差商代替：

（3.2.7）

得到的递推公式：

（3.2.8）

在计算

时，需要用到前两步结果

，称为两步法公式。

（2）积分近似方法

式（3.2.1）的微分方程可写成

，在区间［

上积分，有：

（3.2.9）

上式右边的定积分用不同的积分公方就可得到不同的递推公式。

例如用左矩形公式计算，可以得到显式Euler公式（3.2.4）；用右矩形公式计算，可以得到隐式Euler公式（3.2.6）；

用梯形公式计算．则有递推公式：

（3.2.9）

称为梯形公式。

一般来说，隐式公式的每一次递推计算都需要求解一个非线性方程，虽然可用迭代法求解，但是计算量较大。

为了简化计算过程，可以采用所谓的预测一校正技术，即先用显式公

式计算，得到一个预测值作为隐式公式的迭代初值，然后用隐式公式迭代一次作为非线性方

程的解。

例如梯形公式（3.2.9）可先用显式Euler公式预测，再用梯形公式来校正，即

（3.2.10）

式（3.2.10）称为预测-校正公式或改进Euler公式，也可写成：

（3.2.11）

（3）Taylor展开方法

将初值问题（（3.2.1）在节点

处的函数值

用在节点

处的一阶Taylor展

开式近似表示：

同样可以得到显式Euler公式（3.2.4）。

例题3.1：

数值计算一个粒子在弹性恢复力作用下的位移和速度随时间变化的结果

程序

PROGRAMONE_D_MOTION2

!

Simplestpredictor-corectoralgorithmappliedtoaparticleinone

!

dimensionunderanelasticforce.

USEFGL

USEFML

IMPLICITNONE

INTEGER,PARAMETER:

:

N=101,IN=5

INTEGER:

:

I

REAL:

:

PI,DT

REAL,DIMENSION（N）:

:

T,V,X

PI=4.0*ATAN（1.0）

DT=2.0*PI/100

X

（1）=0.0

T

（1）=0.0

V

（1）=1.0

DOI=1,N-1

T（I+1）=I*DT

!

Predictorforpositionandvelocity

X（I+1）=X（I）+V（I）*DT

V（I+1）=V（I）-X（I）*DT

!

Correctorforpositionandvelocit!

X（I+1）=X（I）+（V（I）+V（I+1））*DT/2.0

V（I+1）=V（I）-（X（I）+X（I+1））*DT/2.0

ENDDO

WRITE（6,"（3F16.8）"）（T（I）,X（I）,V（I）,I=1,N,IN）

CALLMSBACKGROUNDCOLOR（1,1,1）

CALLMSPLOT（MF（T）,MF（X）,'b'）

CALLMSHOLD（'ON'）

CALLMSPLOT（MF（T）,MF（V）,'r'）

CALLMSTITLE（'EULERMETHODEXAMPLE'）

CALLmsYLabel（'X,V'）

CALLmsXLabel（'T'）

CALLMSHOLD（'off'）

CALLMSCAMZOOM（1.2D0）

CALLMSVIEWPAUSE（）

ENDPROGRAMONE_D_MOTION2

Euler方法是常微分方程初值问题数值解法中最简单的一种方法，其精确度较低，但是

可以对它的三种导出方式进行推广，得到更有效的数值方法。

根据递推式中用到的已计算出的数值点个数，求解初值问题的数值方法可以分为单步法

和多步法。

单步法的一般形式可以表示为：

（3.2.12）

即单步法只利用

就可以计算

；如果递推式需要用

的线性组合来计算

，则称为线性多步法，记

，线性多步法可以表示为：

（3.2.13）

单步法只要一个初值

就可以启动递推计算；而式（3.2.13）所示的线性多步法则需要

个初值

，才能够开始递推计算。

根据递推式的右端是否含有待计算的

，求解微分方程的数值方法可以分为显式方法和隐式方法。

如果右端不含

，则只要进行简单递推计算就可得到

，称为显式方法；如果右端也含有

，则必须解方程（组）方可确定

，称为隐式方法。

§3.2.2Runge-Kutta法

Euler方法简单易行，但它的收敛阶数太低。

可以利用Taylor展开式构造高阶的单步方

法。

Taylor展开式实际上是用

在

处的各阶导数值组合来表示

。

所以也属于单步法。

Euler公式可以看成是由一阶Taylor展开式得到的，应用高阶Taylor展开式就可以得到高阶单步法。

例如：

将

在

处作q阶Taylor展开：

（3.2.14）

由于以

满足微分方程，因此它的各阶导数

可以通过函数

对

进行j－1次复合求导获得。

将式（3.2.14）中的余项

舍去，就得到q阶单步方法。

复合函数

的高阶全导数运算，计算量较大，因此高阶Taylor展开方法在实际中难以应用。

Taylor展开式表明，用函数在一个点上的各阶导数值可以近似表示它在邻近另一个点上

的函数值；数值微分公式表明，函数在一个点上的各阶导数值，可以用邻近一些点上的函数

值近似表示。

Runge－Kutta方法，简称RK方法，基于后一种思路，用

在一些特殊点上函数值的线性组合来表示Taylor展开式中的各阶导数值。

N级RK方法的一般形式为：

（3.2.15）

其中：

（3.2.16）

将近似公式和

在

处的Taylor展开式相比较，确定系数，使近似公式具有尽可能高的收敛阶数。

如果式（3.2.15）和式（3.2.16）的局部截断误差:

，则称式（3.2.15）为N级p阶的RK方法。

下面给出2阶RK方法的推导：

取：

（a1）

由

在

的泰勒展开

（a2）

由关系：

，（a2）变为

（a3）

重写（a1）

（a4）

将（a4）的最后一项展开

（a5）