计算方法第八章(常微分方程数值解)
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
下面我们分别取步长为0.1与0.01进行计算， 计算结果显示在下面的图中。
步长为0.1的计算结果。
步长为0.01的计算结果
0.01 0.1 0.2 0.3 0.4 0.41 0.59 0.6 0.9 0.91 0.99 1
0.99005 0.90484 0.81873 0.74082 0.67032 0.66365 0.55433 0.54881 0.40657 0.40252 0.37158 0.36788
我们可得精度更高的欧拉公式： 欧拉中点公式
yn1 yn1 2hf ( xn , yn ) , y( x0 ) y0
y ( xn h ) y ( x n h ) y ( ) 2 y ( xn ) h O (h 2 ) 2h 12
利用中点公式求解微分方程时，有一个问题， 就是计算时需要两个迭代初值！这样的算法 称为二步法。前面的欧拉法称为单步法。 对于这个问题，我们可以先用欧拉公式，通过 给定的初值计算出 的值，然后再利用这两个 值（ y0 和 y1 ）进行计算，直到计算出全部节 点上的值。 一般的单步法： yn1 ( xn , xn1, yn , yn1 ) 一般的 k 步法：
第一节 一般概念 1.1 欧拉法及其简单改进
y( x) f ( x, y ( x)) a x b y (a) y0
方法：选择适当的节点，用差分近似微分，将方程离散化， 从而求在这些节点上的解的近似值。
a x0 x1 x2 xN 1 xN b hn xn1 xn 称为 xn 到 xn1 的步长（通常取为常数h） y ( xn 1 ) y ( xn ) y( x ) |x xn ，记yn为 y(xn）的近似计算值，有 h
欧拉 方法
yn1 yn hf ( xn , yn ) , y( x0 ) y0
例子： y y
y (0) 1
(0 x 1)
f ( x, y ) y yn 1 yn hyn (1 h ) yn y0 1
精确解为 y e x
若用简单迭代，而且只迭代一步，这样组成的一组计算公式称 为预测－－校正公式。（迭代初值 称为校正）
(0) yn 1 yn hf ( xn , yn ) (0) y y [ f ( x , y ) f ( x , y n n n n 1 n 1 )] h / 2 n 1
y( xn ) f ( xn , y( xn )) y ( xn1 ) y ( xn ) hf ( xn , y ( xn )) O(h )
2
而
yn1 y( xn ) hf ( xn , y( xn )) 2 n 1 y ( xn 1 ) yn 1 O(h )
(0) 称为预测，迭代步 yn 1
预测－校正公式也称为改进的欧拉法，将上面的组合公式改 写为：
yn1 yn [ f ( xn , yn ) f ( xn1, yn hf ( xn , yn ))] h / 2
注意到 xn1 xn h ，将上式进一步改写为：
1 yn 1 yn 2 ( K1 K 2 ) K1 hf ( xn , yn ) K hf ( x h, y K ) n n 1 2
例如：后退欧拉法、欧拉梯形公式 显然，利用隐式法求微分方程的数值解是，需要从表达式中 反解未知节点上的函数值。
1.3 隐式法的具体计算： 例如欧拉梯形公式
yn1 yn [ f ( xn , yn ) f ( xn1 , yn1 )] h / 2 yn f ( xn , yn ) h / 2 f ( xn1 , yn1 ) h / 2
作为节点，将被积函数用插值多
项式来近似，用插值多项式带到积分中去求出积分，则可以得
到所谓的亚当斯(Adams)显式公式
h yi 1 yi (b0 fi b1 fi 1 bk fi k ) A
局部截断误差：
R[ y] Bk h
k 1 ( k 1)
y
(i )
类似地，如果取
第八章 常微分方程初值
问题数值解法
本章主要研究常微分方程初值问题的数值求解：
y( x) f ( x, y ( x)) a x b y (a) y0
通常，假设函数 f 关于第二个变量满足李普希茨条件（L条 件），即为存在常数 L > 0，使得
f ( x, y1 ) f ( x, y2 ) L y1 y2
用迭代法计算 yn+1 的值。 （1）简单迭代
(0) yn 1 yn hf ( xn , yn )
( k 1) (k ) yn y f ( x , y ) h / 2 f ( x , y 1 n n n n 1 n 1 ) h / 2
收敛的条件：
Lh / 2 1
用此法解前面的例子
步长0.1
步长0.01
1.4 误差估计 定义：利用第n个节点或之前更多节点的函数精确值，利用近 似公式数值计算第n+1个节点的近似值，所引起的误差，称 为第n+1个节点上的局部截断误差。
我们记 y ( xn 1 ) 为第n+1个节点上解的精确值， yn 1 为假设
yn y ( xn ) 等条件下计算所得的近似值，
xi 1, xi , xi 1,, xi k
作为
节点，可得亚当斯(Adams)隐式公式
h * * * * yi 1 yi (b0 fi 1 b1 fi b2 f i bk f i k ) A
局部截断误差：
R[ y] B h
* k 2 k 1
这是我们最终使用的计算格式。
例子：
y y (0 x 1) y (0) 1
1 yn 1 yn 2 ( K1 K 2 ) K1 hf ( xn , yn ) hyn K hf ( x h, y K ) h( y K ) n n 1 n 1 2
这是所谓的后退欧拉公式。
（3）如选择梯形公式，则得
yn1 yn [ f ( xn , yn ) f ( xn1 , yn1 )] h / 2
这是所谓的欧拉梯形公式。
直接利用已经求得的已知节点上的值计算未知节点上的函数 值的算法称为显式法。 例如：欧拉公式、欧拉中点公式
计算未知节点上的函数值时，用到了未知节点上的函数值， 这种算法称为隐式法。
3
定义：由初值计算到第 n 个节点的近似值与其精确值 之间的误差称为第 n 个节点整体误差。 定理：设下面求解微分方程的数值计算方法
yn1 yn h ( xn , yn , h)
局部截断误差为p+1阶，且函数
( x, y, h) 关于 y 满足
利普希茨条件， ( x, y, h) ( x, y, h) L y y
0.99 0.90438 0.81791 0.7397 0.66897 0.66228 0.55268 0.54716 0.40473 0.40068 0.36973 0.36603
DOUBLE PRECISION h,y(0：100) OPEN(20,FILE='OUTPUT.DAT',STATUS="UNKNOWN") h=1.0/100 y(0)=1.0 do 10 i=1,100 y(i)=y(i-1)*(1.0-h) write(20,*) i*h,y(i) 10 continue END
则局部截断误差为：
n1 y( xn1 ) yn1
p
如局部截断误差为 O(h ) ，称为具有 p 阶局部截断误差。
欧拉方法的误差分析：
y ( xn 1 ) y ( xn ) y ( xn h) y ( xn ) h h y ( xn ) y( xn )h y( )h 2 / 2 y ( xn ) y( xn ) y( )h / 2 h
ynk ( xn ,, xnk , yn ,, ynk 1 , ynk )
1.2 欧拉方法的其他改进 微分方程数值解的关键在于对导数的处理，可以用差分来近似
导数，也可以通过积分，将导数项化掉。
对于方程： y( x) f ( x, y( x)) 首先，作出划分 a x x x x 0 1 2 N 1 xN b 设已经求出第 n 个节点的函数值 yn ，在区间 [ x , x ] 上对 n n 1 方程两边积分
y( xn1 ) y( xn )
积分公式计算积分！
xn1
xn
f ( x, y( x))dx
容易看出，要求第 n+1 个节点的函数值，关键在于选择适当的
（1）如选择下矩形公式，则得
yn1 yn f ( xn , yn )h
这正是前面的欧拉公式。
（2）如选择上矩形公式，则得
yn1 yn f ( xn1 , yn1 )h
同时初值是准确的，则整体截断误差为p阶。 欧拉公式、后退欧拉公式的整体误差为 1 阶。 欧拉中点公式、欧拉梯形公式的整体误差为 2 阶。
微分方程数值解法的进一步改进。再回到恒等式
y( xi 1 ) y( xi )
如果取
xi1
xi
f ( x, y( x))dx
xi , xi 1,, xi k
（2）牛顿迭代
( k 1) yn 1 (k ) (k ) y y f ( x , y ) h / 2 f ( x , y (k ) n 1 n n n n 1 n 1 ) h / 2 yn1 (k ) 1 f y ( xn1 , yn1 ) h / 2
取步长为0.1计算，结果如图。
图：
DOUBLE PRECISION h,y(0:10),ak1,ak2 OPEN(20,FILE='OUTPUT1.DAT',STATUS="UNKNOWN") h=1.0/10 y(0)=1.0 do 10 i=1,10 ak1=-h*y(i-1) ak2=-h*(y(i-1)+ak1) y(i)=y(i-1)+(ak1+ak2)/2.0 10 continue do 20 i=0,10 write(20,*) i*h,y(i),exp(-i*h) 20 continue END
完全类似的可以得到 后退欧拉公式的局部截断误差为：
n1 y( xn1 ) yn1 O(h )
2
欧拉中点公式的局部截断误差为：
源自文库
n1 y( xn1 ) yn1 O(h )
3
欧拉梯形公式的局部截断误差为：
n1 y( xn1 ) yn1 O(h )
同理，对于后退欧拉公式
yn1 yn f ( xn1 , yn1 )h
有预测－校正公式
(0) y n 1 yn hf ( xn , yn ) (0) y y f ( x , y n n 1 n 1 ) h n 1
或改写为：
K hf ( xn , yn ) yn 1 yn f ( xn 1 , yn K )h
y ( xn 1 ) y ( xn ) 欧拉公式中我们利用了近似公式 y ( xn ) h
光这个近似产生的误差为
y ( xn 1 ) y ( xn ) 1 y( xn ) y ( )h O (h ) h 2
y ( xn 1 ) y ( xn 1 ) 利用 y ( xn ) 2h
步长为0.1的计算结果。
步长为0.01的计算结果
0.01 0.1 0.2 0.3 0.4 0.41 0.59 0.6 0.9 0.91 0.99 1
0.99005 0.90484 0.81873 0.74082 0.67032 0.66365 0.55433 0.54881 0.40657 0.40252 0.37158 0.36788
我们可得精度更高的欧拉公式： 欧拉中点公式
yn1 yn1 2hf ( xn , yn ) , y( x0 ) y0
y ( xn h ) y ( x n h ) y ( ) 2 y ( xn ) h O (h 2 ) 2h 12
利用中点公式求解微分方程时，有一个问题， 就是计算时需要两个迭代初值！这样的算法 称为二步法。前面的欧拉法称为单步法。 对于这个问题，我们可以先用欧拉公式，通过 给定的初值计算出 的值，然后再利用这两个 值（ y0 和 y1 ）进行计算，直到计算出全部节 点上的值。 一般的单步法： yn1 ( xn , xn1, yn , yn1 ) 一般的 k 步法：
第一节 一般概念 1.1 欧拉法及其简单改进
y( x) f ( x, y ( x)) a x b y (a) y0
方法：选择适当的节点，用差分近似微分，将方程离散化， 从而求在这些节点上的解的近似值。
a x0 x1 x2 xN 1 xN b hn xn1 xn 称为 xn 到 xn1 的步长（通常取为常数h） y ( xn 1 ) y ( xn ) y( x ) |x xn ，记yn为 y(xn）的近似计算值，有 h
欧拉 方法
yn1 yn hf ( xn , yn ) , y( x0 ) y0
例子： y y
y (0) 1
(0 x 1)
f ( x, y ) y yn 1 yn hyn (1 h ) yn y0 1
精确解为 y e x
若用简单迭代，而且只迭代一步，这样组成的一组计算公式称 为预测－－校正公式。（迭代初值 称为校正）
(0) yn 1 yn hf ( xn , yn ) (0) y y [ f ( x , y ) f ( x , y n n n n 1 n 1 )] h / 2 n 1
y( xn ) f ( xn , y( xn )) y ( xn1 ) y ( xn ) hf ( xn , y ( xn )) O(h )
2
而
yn1 y( xn ) hf ( xn , y( xn )) 2 n 1 y ( xn 1 ) yn 1 O(h )
(0) 称为预测，迭代步 yn 1
预测－校正公式也称为改进的欧拉法，将上面的组合公式改 写为：
yn1 yn [ f ( xn , yn ) f ( xn1, yn hf ( xn , yn ))] h / 2
注意到 xn1 xn h ，将上式进一步改写为：
1 yn 1 yn 2 ( K1 K 2 ) K1 hf ( xn , yn ) K hf ( x h, y K ) n n 1 2
例如：后退欧拉法、欧拉梯形公式 显然，利用隐式法求微分方程的数值解是，需要从表达式中 反解未知节点上的函数值。
1.3 隐式法的具体计算： 例如欧拉梯形公式
yn1 yn [ f ( xn , yn ) f ( xn1 , yn1 )] h / 2 yn f ( xn , yn ) h / 2 f ( xn1 , yn1 ) h / 2
作为节点，将被积函数用插值多
项式来近似，用插值多项式带到积分中去求出积分，则可以得
到所谓的亚当斯(Adams)显式公式
h yi 1 yi (b0 fi b1 fi 1 bk fi k ) A
局部截断误差：
R[ y] Bk h
k 1 ( k 1)
y
(i )
类似地，如果取
第八章 常微分方程初值
问题数值解法
本章主要研究常微分方程初值问题的数值求解：
y( x) f ( x, y ( x)) a x b y (a) y0
通常，假设函数 f 关于第二个变量满足李普希茨条件（L条 件），即为存在常数 L > 0，使得
f ( x, y1 ) f ( x, y2 ) L y1 y2
用迭代法计算 yn+1 的值。 （1）简单迭代
(0) yn 1 yn hf ( xn , yn )
( k 1) (k ) yn y f ( x , y ) h / 2 f ( x , y 1 n n n n 1 n 1 ) h / 2
收敛的条件：
Lh / 2 1
用此法解前面的例子
步长0.1
步长0.01
1.4 误差估计 定义：利用第n个节点或之前更多节点的函数精确值，利用近 似公式数值计算第n+1个节点的近似值，所引起的误差，称 为第n+1个节点上的局部截断误差。
我们记 y ( xn 1 ) 为第n+1个节点上解的精确值， yn 1 为假设
yn y ( xn ) 等条件下计算所得的近似值，
xi 1, xi , xi 1,, xi k
作为
节点，可得亚当斯(Adams)隐式公式
h * * * * yi 1 yi (b0 fi 1 b1 fi b2 f i bk f i k ) A
局部截断误差：
R[ y] B h
* k 2 k 1
这是我们最终使用的计算格式。
例子：
y y (0 x 1) y (0) 1
1 yn 1 yn 2 ( K1 K 2 ) K1 hf ( xn , yn ) hyn K hf ( x h, y K ) h( y K ) n n 1 n 1 2
这是所谓的后退欧拉公式。
（3）如选择梯形公式，则得
yn1 yn [ f ( xn , yn ) f ( xn1 , yn1 )] h / 2
这是所谓的欧拉梯形公式。
直接利用已经求得的已知节点上的值计算未知节点上的函数 值的算法称为显式法。 例如：欧拉公式、欧拉中点公式
计算未知节点上的函数值时，用到了未知节点上的函数值， 这种算法称为隐式法。
3
定义：由初值计算到第 n 个节点的近似值与其精确值 之间的误差称为第 n 个节点整体误差。 定理：设下面求解微分方程的数值计算方法
yn1 yn h ( xn , yn , h)
局部截断误差为p+1阶，且函数
( x, y, h) 关于 y 满足
利普希茨条件， ( x, y, h) ( x, y, h) L y y
0.99 0.90438 0.81791 0.7397 0.66897 0.66228 0.55268 0.54716 0.40473 0.40068 0.36973 0.36603
DOUBLE PRECISION h,y(0：100) OPEN(20,FILE='OUTPUT.DAT',STATUS="UNKNOWN") h=1.0/100 y(0)=1.0 do 10 i=1,100 y(i)=y(i-1)*(1.0-h) write(20,*) i*h,y(i) 10 continue END
则局部截断误差为：
n1 y( xn1 ) yn1
p
如局部截断误差为 O(h ) ，称为具有 p 阶局部截断误差。
欧拉方法的误差分析：
y ( xn 1 ) y ( xn ) y ( xn h) y ( xn ) h h y ( xn ) y( xn )h y( )h 2 / 2 y ( xn ) y( xn ) y( )h / 2 h
ynk ( xn ,, xnk , yn ,, ynk 1 , ynk )
1.2 欧拉方法的其他改进 微分方程数值解的关键在于对导数的处理，可以用差分来近似
导数，也可以通过积分，将导数项化掉。
对于方程： y( x) f ( x, y( x)) 首先，作出划分 a x x x x 0 1 2 N 1 xN b 设已经求出第 n 个节点的函数值 yn ，在区间 [ x , x ] 上对 n n 1 方程两边积分
y( xn1 ) y( xn )
积分公式计算积分！
xn1
xn
f ( x, y( x))dx
容易看出，要求第 n+1 个节点的函数值，关键在于选择适当的
（1）如选择下矩形公式，则得
yn1 yn f ( xn , yn )h
这正是前面的欧拉公式。
（2）如选择上矩形公式，则得
yn1 yn f ( xn1 , yn1 )h
同时初值是准确的，则整体截断误差为p阶。 欧拉公式、后退欧拉公式的整体误差为 1 阶。 欧拉中点公式、欧拉梯形公式的整体误差为 2 阶。
微分方程数值解法的进一步改进。再回到恒等式
y( xi 1 ) y( xi )
如果取
xi1
xi
f ( x, y( x))dx
xi , xi 1,, xi k
（2）牛顿迭代
( k 1) yn 1 (k ) (k ) y y f ( x , y ) h / 2 f ( x , y (k ) n 1 n n n n 1 n 1 ) h / 2 yn1 (k ) 1 f y ( xn1 , yn1 ) h / 2
取步长为0.1计算，结果如图。
图：
DOUBLE PRECISION h,y(0:10),ak1,ak2 OPEN(20,FILE='OUTPUT1.DAT',STATUS="UNKNOWN") h=1.0/10 y(0)=1.0 do 10 i=1,10 ak1=-h*y(i-1) ak2=-h*(y(i-1)+ak1) y(i)=y(i-1)+(ak1+ak2)/2.0 10 continue do 20 i=0,10 write(20,*) i*h,y(i),exp(-i*h) 20 continue END
完全类似的可以得到 后退欧拉公式的局部截断误差为：
n1 y( xn1 ) yn1 O(h )
2
欧拉中点公式的局部截断误差为：
源自文库
n1 y( xn1 ) yn1 O(h )
3
欧拉梯形公式的局部截断误差为：
n1 y( xn1 ) yn1 O(h )
同理，对于后退欧拉公式
yn1 yn f ( xn1 , yn1 )h
有预测－校正公式
(0) y n 1 yn hf ( xn , yn ) (0) y y f ( x , y n n 1 n 1 ) h n 1
或改写为：
K hf ( xn , yn ) yn 1 yn f ( xn 1 , yn K )h
y ( xn 1 ) y ( xn ) 欧拉公式中我们利用了近似公式 y ( xn ) h
光这个近似产生的误差为
y ( xn 1 ) y ( xn ) 1 y( xn ) y ( )h O (h ) h 2
y ( xn 1 ) y ( xn 1 ) 利用 y ( xn ) 2h