大学数值计算方法(第6章 常微分方程数值解法)
常微分方程的数值解法
常微分方程的数值解法在自然科学的许多领域中，都会遇到常微分方程的求解问题。
然而，我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。
在常微分方程课中已经讲过的级数解法，逐步逼近法等就是近似解法。
这些方法可以给出解的近似表达式，通常称为近似解析方法。
还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值。
利用计算机解微分方程主要使用数值方法。
我们考虑一阶常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(yx y y x f dx dy在区间[a, b]上的解，其中f (x, y )为x, y 的已知函数，y 0为给定的初始值，将上述问题的精确解记为y (x )。
数值方法的基本思想是：在解的存在区间上取n + 1个节点b x x x x a n =<<<<= 210这里差i i i x x h -=+1，i = 0,1, …, n 称为由x i 到x i +1的步长。
这些h i 可以不相等，但一般取成相等的，这时na b h -=。
在这些节点上采用离散化方法，（通常用数值积分、微分。
泰勒展开等）将上述初值问题化成关于离散变量的相应问题。
把这个相应问题的解y n 作为y (x n )的近似值。
这样求得的y n 就是上述初值问题在节点x n 上的数值解。
一般说来，不同的离散化导致不同的方法。
§1 欧拉法与改进欧拉法 1.欧拉法1．对常微分方程初始问题(9.2))((9.1) ),(00⎪⎩⎪⎨⎧==y x y y x f dx dy用数值方法求解时，我们总是认为(9.1)、(9.2)的解存在且唯一。
欧拉法是解初值问题的最简单的数值方法。
从（9.2）式由于y (x 0) = y 0已给定，因而可以算出),()('000y x f x y =设x 1 = h 充分小，则近似地有：),()(')()(00001y x f x y hx y x y =≈-(9.3)记 ,n ,,i x y y i i 10 )(== 从而我们可以取),(0001y x hf y y ==作为y (x 1)的近似值。
数值分析常微分方程数值解法
第8页/共105页
➢ 数值积分方法(Euler公式)
设将方程 y=f (x, y)的两端从 xn 到xn+1 求积分, 得
y( xn1) y( xn )
xn1 f ( x, y( x))dx :
xn
xn1 F ( x)dx
xn
用不同的数值积分方法近似上式右端积分, 可以得到计算 y(xn+1)的不同的差分格 式.
h2 2
y''( )
Rn1
:
y( xn1)
yn1
h2 2
y''( )
h2 2
y''( xn ) O(h3 ).
局部截断误差主项
19
第20页/共105页
➢ 向后Euler法的局部截断误差
向后Euler法的计算公式
yn1 yn hf ( xn1, yn1 ), n 0, 1, 2,
定义其局部截断误差为
y 计算 的n递1 推公式，此类计算格式统称为差分格式.
3
第4页/共105页
数值求解一阶常微分方程初值问题
y' f ( x, y), a x b,
y(a)
y0
难点: 如何离散 y ?
➢ 常见离散方法
差商近似导数 数值积分方法 Taylor展开方法
4
第5页/共105页
➢ 差商近似导数(Euler公式)
(0 x 1)
y(0) 1.
解 计算公式为
yn1
yn
hfn
yn
h( yn
2xn ), yn
y0 1.0
n 0, 1, 2,
取步长h=0.1, 计算结果见下表
13
常微分方程的数值解法
第六章 常微分方程的数值解法 §6.0 引言§6.1 算法构造的主要途径§6.2 Runge-Kutta Method 算法§6.3 线性多步法§6.4 线性多步法的一般形式§6.5 算法的稳定性、收敛性§6.0 引 言1． 主要考虑如下的一阶常微分方程初值问题的求解:()()00,dy f x y dx y x y ⎧=⎪⎪⎨=⎪⎪⎩ 微分方程的解就是求一个函数y=y(x)，使得该函数满足微分方程并且符合初值条件。
2. 例如微分方程:xy '-2y=4x ；初始条件: y(1)=-3。
于是可得一阶常微分方程的初始问题24(1)3y y x y ⎧'=+⎪⎨⎪=-⎩。
显然函数y(x)=x 2-4x 满足以上条件,因而是该初始问题的微分方程的解。
3. 但是,只有一些特殊类型的微分方程问题能够得到用解析表达式表示的函数解，而大量的微分方程问题很难得到其解析解，有的甚至无法用解析表达式来表示。
因此，只能依赖于数值方法去获得微分方程的数值解。
4． 微分方程的数值解：设微分方程问题的解y(x)的存在区间是[a,b ]，初始点x 0=a ，将[a,b ]进行划分得一系列节点x 0 , x 1 ,...,x n ，其中a= x 0< x 1<…< x n =b 。
y(x)的解析表达式不容易得到或根本无法得到，我们用数值方法求得y(x)在每个节点x k 的近似值y(x k )，即y≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。
如图所示:§6.1 算法构造的主要途径x 0 x 1 x 2 ...1 欧拉公式1.1 构造的思想：利用差商代替一阶导数，即010()()x x y x y x dy dx h =-≈，则 1000()()(,)y x y x f x y h -≈。
常微分方程数值解法
介绍常微分方程数值解法常微分方程（ordinary differential equations，ODE）可用于描述许多日常存在的物理系统。
处理ODE问题常常被称为数值求解法，这指的是找到概括ODE或者其他适用于数学模型的解决方案来模括这些ODE。
这种解决方案可能在一系列不同方案中发挥重要作用，以此来提供更好的解释和预测。
常微分方程与几何图形更为相关，它利用二维或者三维空间中曲线的绘制以及分析。
通过引入一些不同的方法，可以对不同的常微分方程中的量进行描述，使得可以通过数值方法的解析来进行研究。
数值解法可能是时间消耗较多的，但有助于验证几何图形中的某些过程，以此帮助揭示数学模型。
四种常见的常微分方程数值解法四种常见的常微分方程数值解法是：前向差分法、向后差分法、中点法和全分方法。
•前向差分法：前向差分法的基本概念是利用ODE的特定解来表达时间步的影响。
这是一种基本的数值法，可以在ODE中确定任意位置的点作为终点。
在这里，任何这样的点都可以表示为ODE右边的时间步。
•向后差分法：它是反过来基于前向差分法。
它要求对ODE中的时间步进行逆向推导，以获得某一特定点的解。
向后差分法要求推导反向解中点，以便可以从每一步中获取该点的解。
•中点法：这是一种非常基本的数值解法，可以用来求解ODE中的某一步的解，但不具有直观的方法解释。
主要的思想是在每一次时间步中通过求出ODE的中点来寻找解。
•全分方法：这是一种更复杂的数值解法，它要求将ODE中的每一步解细分并解决。
与前面提到的三种解法不同，它首先要求将ODE分解成若干离散区间，然后计算每一段区间中的点。
这种解法可以更准确地进行处理，但时间消耗较多，因此比较少被使用。
优化方案在需要解决常微分方程时，为了得到最佳的结果，有必要考虑一些优化措施。
•首先，应考虑将一个复杂的ODE拆分成一些更易解决的问题。
这样做的结果是，预见到解决此ODR的总耗时将会降低。
•其次，为了加快计算速度，可以考虑使用预解算法。
哈尔滨工程大学工程算法课件06常微分方程的数值求解
欧拉法得： yn 1 yn hf xn , yn 因此，局部截断误差是 o h 2 。
19
2 改进Euler法
2.1方法构造
dy f x, y ，对其从 xk 到 xk 1 进 在微分方程初值问题 dx 行定积分得：
y xk 1 y xk
yk 1 是未知,待求的,未知量在 f x, y 中这是
一个方程,如f是非线性或超越函数,此方程是无法直接解出来(要 依靠迭代法才能解出)。这类格式称为隐式格式。
21
2.3 算例
y y x 例：用改进欧拉公式求解 , h 0.2 y 0 2 解： f x, y y x h yk 1 yk f xk , yk f xk 1 , yk 1 2 h h 1 2 y 2 x x y k 1 k 1 h k h k 1 1 2 2 可以从隐式格式中解出 yk 1 问题的精确解是 y x e x x 1
16
精确解为： y x 2 x
2
可以看出误差随着计算在积累。
17
1.4 Euler法的特点和误差
迭代格式 特点
1 单步方法：
yn 1 yn hf xn , yn n 0,1, 2,, N 1
2 显示格式： 3 局部截断误差为O h2
18
第六章 常微分方程数值解
§6.0 引言
§6.1 欧拉方法 §6.2 龙格-库塔方法
§6.3 单步法的收敛性和稳定性
§6.4 线性多步法
1
§6.0 引言
1 主要考虑如下的一阶常微分方程初值问题 的求解：
dy f x, y dx y x0 y0
实验09 数值微积分与方程数值解(第6章)
实验09 数值微积分与方程数值求解（第6章 MATLAB 数值计算）一、实验目的二、实验内容1. 求函数在指定点的数值导数232()123,1,2,3026x x x f x x xx x==2. 用数值方法求定积分(1) 210I π=⎰的近似值。
程序及运行结果：《数学软件》课内实验王平(2) 2221I dx x π=+⎰程序及运行结果：3. 分别用3种不同的数值方法解线性方程组6525494133422139211x y z u x y z u x y z u x y u +-+=-⎧⎪-+-=⎪⎨++-=⎪⎪-+=⎩ 程序及运行结果：4. 求非齐次线性方程组的通解1234123412342736352249472x x x x x x x x x x x x +++=⎧⎪+++=⎨⎪+++=⎩5. 求代数方程的数值解(1) 3x +sin x -e x =0在x 0=1.5附近的根。
程序及运行结果（提示：要用教材中的函数程序line_solution ）：(2) 在给定的初值x 0=1，y 0=1，z 0=1下，求方程组的数值解。
23sin ln 70321050y x y z x z x y z ⎧++-=⎪+-+=⎨⎪++-=⎩6. 求函数在指定区间的极值(1) 3cos log ()xx x x xf x e ++=在(0,1)内的最小值。
(2) 33212112122(,)2410f x x x x x x x x =+-+在[0,0]附近的最小值点和最小值。
7. 求微分方程的数值解，并绘制解的曲线2250(0)0'(0)0xd y dyy dx dx y y ⎧-+=⎪⎪⎪=⎨⎪=⎪⎪⎩程序及运行结果（注意：参数中不能取0，用足够小的正数代替）：令y 2=y,y 1=y '，将二阶方程转化为一阶方程组：'112'211251(0)0,(0)0y y y x x y y y y ⎧=-⎪⎪=⎨⎪==⎪⎩8. 求微分方程组的数值解，并绘制解的曲线123213312123'''0.51(0)0,(0)1,(0)1y y y y y y y y y y y y =⎧⎪=-⎪⎨=-⎪⎪===⎩程序及运行结果：三、实验提示四、教程：第6章 MATLAB 数值计算（2/2）6.2 数值微积分 p155 6.2.1 数值微分1. 数值差分与差商对任意函数f(x)，假设h>0。
常微分方程的数值解法
常微分方程的数值解法常微分方程是研究变量的变化率与其当前状态之间的关系的数学分支。
它在物理、工程、经济等领域有着广泛的应用。
解常微分方程的精确解往往十分困难甚至不可得，因此数值解法在实际问题中起到了重要的作用。
本文将介绍常见的常微分方程的数值解法，并比较其优缺点。
1. 欧拉方法欧拉方法是最简单的数值解法之一。
它基于近似替代的思想，将微分方程中的导数用差商近似表示。
具体步骤如下：（1）确定初始条件，即问题的初值。
（2）选择相应的步长h。
（3）根据微分方程的定义使用近似来计算下一个点的值。
欧拉方法的计算简单，但是由于误差累积，精度较低。
2. 改进欧拉方法为了提高欧拉方法的精度，改进欧拉方法应运而生。
改进欧拉方法通过使用两个点的斜率的平均值来计算下一个点的值。
具体步骤如下：（1）确定初始条件，即问题的初值。
（2）选择相应的步长h。
（3）根据微分方程的定义使用近似来计算下一个点的值。
改进欧拉方法相较于欧拉方法而言，精度更高。
3. 龙格-库塔法龙格-库塔法（Runge-Kutta）是常微分方程数值解法中最常用的方法之一。
它通过迭代逼近精确解，并在每一步中计算出多个斜率的加权平均值。
具体步骤如下：（1）确定初始条件，即问题的初值。
（2）选择相应的步长h。
（3）计算各阶导数的导数值。
（4）根据权重系数计算下一个点的值。
与欧拉方法和改进欧拉方法相比，龙格-库塔法的精度更高，但计算量也更大。
4. 亚当斯法亚当斯法（Adams）是一种多步法，它利用之前的解来近似下一个点的值。
具体步骤如下：（1）确定初始条件，即问题的初值。
（2）选择相应的步长h。
（3）通过隐式或显式的方式计算下一个点的值。
亚当斯法可以提高精度，并且比龙格-库塔法更加高效。
5. 多步法和多级法除了亚当斯法，还有其他的多步法和多级法可以用于解常微分方程。
多步法通过利用多个点的值来逼近解，从而提高精度。
而多级法则将步长进行分割，分别计算每个子问题的解，再进行组合得到整体解。
第6章常微分方程初值问题的解法
ykh 2 k[ (ykx k 1 ) ( yk 1x k 1 1 )]
yk11 29 1yk1k05110
预估－校正Euler方法：
y k 1 0 .90 y k 5 0 .00 k 9 0 .1 5
20
Euler方法
xk
yk
yk y(xk)
0.0 1.000000
0.0
梯形方法
yk
yk y(xk)
1.000000
0.0
续
预估－校正方法
yk
yk y(xk)
1.000000
0.0
0.1 1.000000 0.2 1.010000
4.8×10－3 8.7×10－3
1.004762 1.018594
y(0) 1
其解析解为： y1xe-t2dt x[0,1] 0 很难得到其解析解
4
例如：
y=x+y , x[0,1]
y(0) 1
其解析解为 yx12ex
只有一些特殊类型的微分方程问题能够得到用解析表达式 表示的函数解，而大量的微分方程问题很难得到其解析解。
因此，只能依赖于数值方法去获得微分方程的数值解。
例如：
y=x+y , x[0,1]
y(0) 1
其解析解为：yx12ex
3
但是, 只有一些特殊类型的微分方程问题能够得到用解析 表达式表示的函数解，而大量的微分方程问题很难得到其解 析解。
因此，只能依赖于数值方法去获得微分方程的数值解。
例如：
y =e-x2 ,
x[0,1]
7.5×10－5 1.4×10－4
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
6.1 初值问题的Euler方法
设一阶常微分方程初值问题 dy f ( x, y ) dx y ( x0 ) y0 记 xn a nh (n 0,1,2,...), h 0为步长，一般总 假定为常数。该式的数值解是指通过某种方法 去获得解y ( x)在点xn 上的近似值yn , 即 y ( xn ) y n (n 0,1,2,...)
误差概述
一般说来，一个算法，局部截断 误差阶ｐ 越大， 则精度相对的越高。 h2 对于显式Euler法：Tn 1 y( xn ) O(h 3 ) 2 h2 对于隐式Euler法：Tn 1 y( xn ) O( h 3 ) 2 h3 ( xn ) O (h 4 ) 对于梯形公式：Tn 1 y 2
这就是隐式的`Euler 公式或向后Euler方法，它与显式 的不同在于，它每算一步要解函数方程(2)才能得到 yn 1。
初值问题的Euler方法
如果取以上两式的算术平均值的结果，则得 h yn 1 yn [ f ( xn , yn ) f ( xn 1 , yn 1 )] （n 0,1,2,...） 2 称为梯形公式。 计算yn时常用以下迭代式：
~ ~
误差概述
定义6.1. 如果给果给定的算法的截断 2 误差为 Tn 1 O(h p 1 )
p 1
则称该算法具有p阶精度。如果 Tn 1 g ( xn ,y ( xn )) h O(h
p2
)
则非零项g ( xn , y ( xn )) h p 1称为为局部截断误差主。
1 yn 1 yn (k1 k 2 ) 2 2 xn 改进的Euler法： k1 0.1( yn ) yn 2( xn 0.1) k 2 0.1( yn k1 ) yn k1 计算结果如下表所示：
(n 0,1,2,...)
初值问题的Euler方法
对于( 2)计算yn 1 ,由于迭代工作量较大, 一般只 迭代一次, 构成一类预估 校正算法, 即
( p) yn 1 yn hf ( xn , yn ) h (c) ( p) yn ` yn [ f ( xn , yn ) f ( xn 1 , yn 1 )] 2 ( 并取yn 1 ync )1。
6.1.2 误差概述
显式单步法一般形式为 yn 1 yn h ( xn , yn , h) 而隐式单步法一般形式为 yn 1 yn h ( xn , yn , xn 1 , yn 1 , h) 函数与f ( x, y )有关，称为增量函数。
误差概述
定义6.1.1 从初值y ( x0 ) y0出发，由单步法显式
或隐式逐步计算，得xn 1的值yn 1 , 则en 1 y ( xn 1 ) yn 1 称为在点xn 1上的整体截断误差。如果第n步在点xn的 值计算没有误差，即yn y ( xn ),由单步法计算出 y n 1 , 则Tn 1 y ( xn 1 ) y n 1 , 称为点xn 1上的局部截断误差。
于是 yn 1 y ( xn ) (1 2 )hy( xn ) [ 2 f x( xn , yn ) 2 f ( xn , yn ) f y ( xn , yn )]h 2 O(h3 )
Runge-Kutta方法
与Taylor展式相比较得 : 1 2 1 1 2 2 2 1 2 由于四个参数，三个方程，因此有一个 自由参数，即解答不唯一。
这就是显式的`Euler公式，它可以从y0出发，逐次
初值问题的Euler方法
如果用xn 1代替x0 , 于是该式可离散为： y ( xn h ) y ( xn ) f ( xn 1 , y ( xn 1 )) h 以yn 表示y ( xn )的近似值，则有 yn 1 yn hf ( xn 1 , yn 1 ) （n 0,1,2,...） ( 2)
初值问题的Euler方法
例6.1.1 设初值问题
2x dy y y dx y (0) 1 试分别用Euler法和改进Euler法求解, 并与 精确解y 1 2 x 进行比较。
初值问题的Euler方法
解：取h 0.1, 计算x [0,1]上结果，此时 2 xn Euler法：yn 1 yn 0.1( yn ) yn (n 0,1,2,...)
( yn0 )1 yn hf ( xn , yn ) （3） h ( k 1) (k ) yn 1 yn 2 [ f ( xn , yn ) f ( xn 1 , yn 1 )] （k 0,1,2,...） ( ( 当 | ynk 1) ynk ) | 1 ( 时, 取yn 1 ynk 1) 1
证明：由式（2）和（3）有 | y n 1 y
(k 1) n 1
h ( ) | | f ( xn 1 , yn 1 ) f ( xn 1 , ynk 1 ) | 2 hL ( ) | yn 1 ynk 1 | 2 ......
hL k 1 ( ( ) | yn 1 yn0 )1 | 2 hL k 1 ( 由假设知 : lim ( ) 0 , 故有 lim ynk 1) yn 1。 1 k 2 k
初值问题的Euler方法
定理6.1.1 设函数f ( x, y )对变量y满足Lipschitz
hL 条件，L为Lipschitz常数。如果步长h满足0 1, 2 2 ( ) 即h 时，则由（3）产生的序列{ ynk 1}( k 0,1,2...) L 收敛。
初值问题的Euler方法
1 1 (1) 取1 , 可得2 , 1, 此时算式为 2 2 1 yn 1 yn 2 (k1 k 2 ) k1 hf ( xn , yn ) k hf ( x h, y k ) n n 1 2 这是改进的Euler方法。
~ ~
数值稳定性分析
定义6.1.3 若某数值算法的绝对稳定性区 域包含hλ平面上的左半平面Re(hλ)＜0， 则称该方法是A稳定的。 隐式Euler法是A稳定的。
6.2 Runge-Kutta方法
受改进的Euler方法启发，更一般算式可设为 yn 1 yn 1k1 2 k 2 (n 0,1,2,...) k1 hf ( xn , yn ) k hf ( x h, y k ) n n 1 2 适当选择参数1，2，，，使局部截断误差 Tn 1 y ( xn 1 ) yn 1 O(h ), 这里仍假定yn y ( xn )。
1 3 2 (3) 取1 , 可得2 , , 又有算式 4 4 3 1 （ yn 1 yn 4 k1 3k 2 ) k1 hf ( xn , yn ) 2 2 k 2 hf ( xn h, yn k1 ) 3 3 这也是二阶R - K方法。
初值问题的Euler方法
上式还常写成 1 yn 1 yn 2 f (k1 k 2 ) k1 hf ( xn , yn ) k hf ( x h, y k ) (n 0,1,2,...) n n 1 2 该式称为改进Euler方法, 亦可写成 h yn 1 yn [ f ( xn yn ) f ( xn 1 , yn hf ( xn , yn ))] 2
Runge-Kutta方法
1 (2) 取1 0, 可得2 1, , 此时算式为 2 yn 1 yn k 2 k1 hf ( xn , yn ) 1 1 k 2 hf ( xn h, yn k1 ) 2 2 这是二阶R - K方法.
初值问题的Euler方法
为实现这一目标，Euler方法首先将微分算子离 散化，并用xn 代替x0 , 于是该式可离散为： y ( xn h ) y ( xn ) f ( xn , y ( xn )) h 以yn 表示y ( xn )的近似值，则有 yn 1 yn hf ( xn , yn ) （n 0,1,2,...） 算出y1 , y2 , y3 ...。 (1)
3
Runge-Kutta方法
由二元函数Taylor展开式 : k 2 hf ( xn , yn ) h 2 f x( xn , yn ) hk1 f y ( xn , yn ) O(h 3 ) hy( xn ) h 2 (f x( xn , yn ) f ( xn , yn ) f y ( xn , yn )) O(h 3 )
x
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.000 1.000000 1.191818 1.277438 1.358213 1.435133 1.508966 1.580338 1.649783 1.717779 1.784770
改进的Euler法y 精确解
1.000000 1.095909 1.184097 1.266201 1.343360 1.416402 1.485956 1.552514 1.616475 1.678166 1.737867 1.000000 1.095445 1.183216 1.264911 1.341641 1.414214 1.483240 1.549193 1.612452 1.673320 1.732051
第6章 常微分方程数值解法
绪论
在工程和科学计算中，所建立的各 种常微分方程的初值或边值问题，除很 少几类的特殊方程能给出解析解，绝大 多数的方程是很难甚至不可能给出解析 解的，其主要原因在于积分工具的局限 性。因此，人们转向用数值方法去解常 微分方程，并获得相当大的成功，讨论 和研究常微分方程的数值解法是有重要 意义的。