常微分方程数值解法

第九章 常微分方程的数值解法主要内容§1、引言

§2、初值问题的数值解法--单步法§3、龙格-库塔方法

§4、收敛性与稳定性§5、初值问题的数值解法―多步法 §6、方程组和刚性方程 §7、习题和总结

§1、 引

言

主要内容

研究的问题 数值解法的意义

1.什么是微分方程 ? 现实世界中大多数事物 内部联系非常复杂

其状态随着 时间、地点、条件 的不同而不同 找出其状态和状态变化规律之间的相互联系， 也即一个或一些函数与他们的导数之间的关系

此种关系的数学表达就为

微分方程

2.数值求解微分方程的意义如何建立数学模型已在建模课程中得到讨论， 各类微分方程本身和他们的解所具有的特性 已在常微分方程及数学物理方程中得以解释， 本章专门 讨论

如何利用数值方法求解微分方程(组)的问题。

3.什么是微分方程 (组)的解析解?

3.什么是微分方 程(组)的解析解?

一个具有所要求阶连续导数的解析函数，将 它代入微分方程，恰使其所有条件都得到满 足的解称为解析解(或古典解),称为真解或解。 寻找解析解的过程称为求解微分方程。 y f ( x , y ), y ( x 0 ) y0 . (1.1) (1.2)

虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问 题中归结出来的微分方程主要靠数值解法. 所谓数值解法, 就是寻求解y(x)在一系列离散节点

x1 x2 xn xn 1 上的近似值 y1,y2, ,yn,yn+1, . 相邻两个节点的间距

hn=xn+1-xn称为步长. 今后如不特别说明，总是假定hi=h(i=1,2, )为定数, 这时节点为xn=x0+nh(i=0,1,2, ) (等距节点).

初值问题的数值解法有个基本特点，他们都采 取“步进式”,即求解过程顺着节点排列的次序一步 一步地向前推进. 描述这类算法，只要给出用已知信

息yn,yn-1,yn-2, 计算yn+1的递推公式.首先，要对微分方程离散化，建立求解数值解的

递推公式. 一类是计算yn+1时只用到前一点的值yn,称为单步法. 另一类是用到yn+1前面 k 点的值yn,yn-1, ,

yn-k+1,称为k步法. 其次，要研究公式的局部截断误差和阶，数值解yn与精确解y(xn)的误差估计及收敛性， 还有递推公式的计算稳定性等问题.

5.常微分方程数值解法的特点常微分方程的数值解法常用来求近似解 根据提供的算法 通过计算机

数值解法得到的近似解(含误差)是一个 离散的函数表.

便捷地实现

6.基本知识本章主要讨论一阶常微方程的初值问题 各 种 数 值 解 法

y f ( x, y ) y( x0 ) y0

a x b

(1.1) (1.2)

其中f (x,y)是已知函数，(1.2)是定解条件也称为 初值条件。

§2、初值问题的数值解法―

单步法

简单的数值方法与基本概念1. 简单欧拉法(Euler)

2.后退的欧拉法3.梯形法 4.改进Euler法

1. 简单的欧拉(Euler)方法考虑模型：

f ( x, y) y y( x0 ) y0欧 拉 方 法

a x b

(1.1) (1.2)

最简单而直观 实用方法弄清常微方程初值 问题数值解法的一 些基本概念和构造 方法的思路.

在精度要求不高时 通过欧拉方法的讨论

2. 欧拉方法的导出把区间[a,b] 分为n个小区间N等分

步长为 hi ( xi 1 - xi )节点 xi a ihi , 一般取hi h( (b a) / n)即等距

要计算出解函数 y(x) 在一系列节点a x0 x1 xn b

处的近似值

yi y( xi )

y f ( x, y ) y ( x0 ) y0

a x b

(1.1) (1.2)

对微分方程(1.1)两端从 xn到xn 1 进行积分

xn 1

xn

y dx

xn 1

xn

f ( x, y( x))dxxn 1

y( xn 1 ) y( xn )

xn

f ( x, y( x))dx

令g ( x) f ( x, y( x))

右端积分用 左矩形数值 求积公式

b

a

g ( ) g ( x) dx (b a ) g ( a ) (b a ) 2 2

得

yn 1 yn ( x

x ) f ( xn , y( xn )) n 1 n hf ( x , y( x )) n n

或用向前差 商近似导数

y ( xn 1 ) y ( xn ) y ( xn ) hx0 x1

y( xn 1 ) y( xn ) hy ( xn ) yn h f ( xn , yn )

yi 1 yi h f ( xi , yi ) (i 0, ... , n 1)

依上述公式逐次计算可得：y1 y0 hf ( x0 , y0 ) y2 y1 hf ( x1 , y1 ) yn 1 yn hf ( xn , yn )亦称为欧拉折线法/* Euler’s polygonal arc method*/

每步计算 yn 1只用到 yn 故也称Euler为单步法。公式右端只含有已知项 所以又称为显格式的单步法。

例1 用欧拉公式求解初值问题 2x (0 x 1), y y y y(0) 1.

( 2.2)

解 取步长h=0.1,欧拉公式的具体形式为 2 xn yn 1 yn h( yn ) yn 其中xn=nh=0.1n (n=0,1, ,10), 已知y0 =1, 由此式可得 2 x0 y1 y0 h( y0 ) 1 0.1 1.1 y0

2 x1 0.2 y2 y1 h( y1 ) 1.1 0.1(1.1 ) 1.191818 y1 1.1

依次计算下去，部分计算结果见下表.xn 欧拉公式数值解yn 准确解y(xn) 0.2 0.4 0.6 0.8 1.0 1.191818 1.358213 1.508966 1.649783 1.784770 1.183216 1.341641 1.483240 1.612452 1.732051 误差 0.008602 0.016572 0.025726 0.037331 0.052719

与准确解 y 1 2 x 相比，可看出欧拉公式的计算结 果精度很差.

欧拉公式具有明显的几何意义, 就是用折线近似 代替方程的解曲线，因而常称公式(2.1)为欧拉折线法. 还可以通过几何直观来考察欧拉方法的精度.假 设yn=y(xn),即顶点Pn落在积分曲线y=y(x)上，那么， 按欧拉方法做出的折线 PnPn+1便是y=y(x)过点Pn 的切线.从图形上看,这 样定出的顶点Pn+1显著 地偏离了原来的积分曲 线，可见欧拉方法是相 当粗糙

的.

ypn

pn 1y y ( x)p xn 1

xn

xn 1

x

定义 在假设 yi = y(xi),即第 i 步计算是精确的前提下，考 虑的截断误差 Ri = y(xi+1) yi+1 称为局部截断误差 /* local truncation error */。 定义 若某算法的局部截断误差为 O(hp+1),则称该算法有 p 阶精度。

欧拉法的局部截断误差：Ri y( xi 1 ) yi 1 [ y( xi ) hy ( xi ) h2 y ( xi ) O(h3 )] [ yi hf ( xi , yi )]2

Ri 的主项

/* leading term */

h2 y ( xi ) O(h3 )2

欧拉法具有 1 阶精度。

5. 欧拉公式的改进：

隐式欧拉法 /* implicit Euler method */y ( x1 )

向后差商近似导数

y ( x1 ) y0 h f ( x1 , y ( x1 ))

y( x1 ) y( x0 ) h

x0

x1

y i 1 y i h f ( x i 1 , yi 1 ) ( i 0, ... , n 1)

由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故 称为隐式 /* implicit */ 欧拉公式，而前者称为显式 /* explicit */ 欧拉公式。 一般先用显式计算一个初值，再迭代求解。

隐式欧拉法的局部截断误差：Ri y( xi 1 ) yi 1 h2 y ( xi ) O(h3 )2

即隐式欧拉公式具有 1 阶精度。

设用欧拉公式

y

(0) n 1

yn hf ( xn , yn )

给出迭代初值 y

(0) n 1 ,用它代入(2.5)式的右端，使之转

化为显式，直接计算得

y然后再用 y

(1) n 1代入(2.5)式，又有 ( 2) (1) n 1 n n 1 n 1

(1) n 1

yn hf ( xn 1 , y ),(0) n 1

y

y hf ( x , y ).(k ) n 1

如此反复进行，得

y

( k 1) n 1

yn hf ( xn 1 , y ) (k 0,1, ). (2.6)

本文来源：https://www.bwwdw.com/article/ytwe.html