matlab软件欧拉算法教程
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
改进的Euler方法
第二节 改进的Euler方法
已知初值问题的一般形 式为: dy f ( x, y) dx y ( x 0 ) y0 a xb (1)
将方程两端从xn到xn+1求积分，得 y( xn1 ) y( xn )
xn1
xn
f ( x , y( x ))dx
预测—校正系统称作改进的欧拉公式。
单步显式格式
改进的Euler方法
改进Euler格式的嵌套形式
h yn1 yn f ( xn , yn )+f ( xn1 , yn +hf ( xn , yn ) 2
平均化形式： y p yn hf ( xn , yn ) yc yn hf ( xn 1 , y p ) 1 yn 1 ( y p yc ) 2 y ( x 0 ) y0
Step 1: 将 K2 在 ( xn , yn ) 点作 Taylor 展开
K 2 f ( xn ph , yn phK1 ) f ( xn , yn ) phf x ( xn , yn ) phK1 f y ( xn , yn ) O( h2 )
y( xn ) phy( xn ) O(h2 )
四、四阶Runge-Kutta方法
继续上述过程，可以进一步导出四阶Runge-Kutta 格式
yn1 yn h K1 2 K 2 2 K 3 K 4 , 6 K f x , y , n n 1 h K 2 f x n 1 , yn K1 , 2 2 K f x , y h K , n 1 n 2 3 2 2 K 4 f xn1 , yn hK 3 .
经典R-K公式 每一步计算需 要四个函数值
注意的问题
Runge-Kutta法的主要运算在于计算 Ki 的值，即计算 f 的值。计算量与可达到的最高精度阶数的关系：
每步须算Ki 的个数
2 3 4 5 6 7
n8
O(hn2 )
4 3 4 5 可达到的最高精度 O(h2 ) O(h ) O(h ) O(h ) O(h )
等于中点的斜率值f xn 1 , y( xn 1 ) , 本格式表示用中
2 2
点的斜率代替平均斜率K*
计算量：每步需要计算两次函数f的值
Q: 为获得更高的精度，应该如何进一步推广？
三、三阶Runge-Kutta方法
为进一步提高精度，设除xn+p外再考察一点
xnq xn qh,
p q 1
利用数值积分求积分项的方法离散
(1)左矩形法
xn1
ຫໍສະໝຸດ Baidu
f ( x , y( x ))dx hf ( xn , y( xn )) O( h2 )
xn
y( xn1 ) y( xn ) hf ( xn , y( xn ))
Euler格式 一阶方法
改进的Euler方法
(2)梯形法
xn1
O(h6 )
R-K(高阶)方法不唯一,选择不同的参数能得到 不同的R-K公式 R-K方法的推导是基于Taylor展开法，因而要求 解具有较好的光滑性，如果光滑性较差精度可 能不如改进Euler方法,最好采用低阶算法而将 步长h 取小。
四阶R-K方法实现
开始
输入 x 0, y 0, h, N
n N 1; x1 x0 , y1 y0
xn
h f ( x , y( x ))dx f ( xn , y( xn )) f ( xn1 , y( xn1 )) O( h3 ) 2
h yn1 yn [ f ( xn , yn ) f ( xn1 , yn1 )] 梯形格式 2
梯形格式是显式Euler格式与隐式Euler格式的 算术平均 Euler格式是显式算法，计算量小，但精度低 梯形格式，精度较高，但是隐式算法，需要通过 迭代过程求解，计算量大
斜率 一定取K1 K2 的平均值吗？
步长一定是一个h 吗？
用xn , xn1两个点的斜率值K1 , K 2的算术平均作为平均 斜率K *；xn 1处的斜率值K 2，由已知信息y n 通过Euler 公式预报
Runge-Kutta方法的设计思想 设法在[xn,xn+1]区间内多预报几个点的斜率值， 利用这些斜率值，将他们加权平均作为平均斜率 的近似，有可能构造出更高精度的计算格式
y( xn1 ) y( xn ) hK
*
寻求计算平均斜率的算法
考察欧拉法，以xn的斜率值
K1 f ( xn , yn )
作为平均斜率
考察改进的欧拉法，可以将其改写为：
1 1 yi 1 yn h K1 K 2 2 2 K1 f ( x n , yn ) K 2 f ( xn h, yn hK1 )
ynq yn qh(1 )K1 K2
K 3 f ( x n q , yn q )
因此K 3为:
利用Taylor展开法选择参数p, q, , , ，使此计算格式 具有三阶精度这类格式统称为三阶Runge-Kutta格式
R-K法的常用公式
常用的三阶R-K方法.
3 要求 Rn y( xn1 ) yn1 O(h ) ，则必须有：
1 1 2 1 , 2 p 2
这里有 3 个未知 数， 2 个方程。
存在无穷多个解。所有满足上式的格式统称为2阶 龙格 - 库塔格式。
1 (1) p 1, 1 2 为改进Euler 格式 2 1 (2)p , 1 0, 2 1为变形的Euler 格式 中点格式 2 yn h K1是Euler 法预报的中点xn 1 的近似值，K 2 近似 2 2
h yn 1 yn 6 K1 4 K 2 K 3 , K1 f x n , yn , K 2 f x 1 , yn h K1 , 2 n+ 2 K 3 f xn +1 , yn +h( K1 2 K 2 ) .
二、二阶Runge-Kutta方法
推广改进的Euler方法，在区间[xn , xn1 ]内任取一点 xn p xn ph,
* *
0 p 1
希望用xn，xn p两个点的斜率K1 , K 2加权平均得到平 均斜率K ，即 K =1 K1 2 K 2 , 1 ,2为待定参数
Y
x1 b
N 结束
例题3
例2
求解初值问题 步长h 0.1) ( 2x y y y y ( 0) 1 (0 x 1)
解
f ( x, y) y 2 x / y
y p yn hf ( xn , yn ) 改进Euler格式 yc yn hf ( xn1 , y p ) yn1 1 ( y p yc ) 2
Step 3: 将 yn+1 与 y( xn+1 ) 在 xn 点的泰勒展开作比较
yn1 yn (1 2 )h y( xn ) 2 ph2 y( xn ) O(h3 )
h2 y( xn1 ) y( xn ) hy( xn ) y( xn ) O( h3 ) 2
Step 2: 将 K2 代入第1式，得到
yn1 yn h 1 y( xn ) 2 [ y( xn ) phy( xn ) O( h2 )] yn (1 2 )h y( xn ) 2 ph2 y( xn ) O( h3 )
§2 Runge-Kutta Method
Euler法
xn Yn
0.1 1.1 0.2 1.1918
改进Euler法
0.0046 1.0959 0.0086 1.1841 0.0005 0.0009
准确解
1.0954 1.1832
|yn-y(xn)| Yn
|yn-y(xn)| y(xn)
0.3 1.2774
0.4 1.3582 0.5 1.4351 0.6 1.5090 0.7 1.5803 0.8 1.6498 0.9 1.7178 1.0 1.7848
改进Euler方法计算框图 开始
输入x0 , y0 , h, b
n 1
x1 x0 h y p y0 hf ( x0 , y0 ) yc y0 hf ( x1 , y p ) 1 y1 ( yc y p ) 2
输出x1 , y1
n n 1 x0 x1 , y0 y1
0.0125 1.2662
0.0166 1.3434 0.0209 1.4164 0.0257 1.4860 0.0311 1.5525 0.0373 1.6165 0.0445 1.6782 0.0527 1.7379
0.0013
0.0018 0.0022 0.0028 0.0033 0.0040 0.0049 0.0058
用三个点x n , x n p , x nq的斜率K 1 , K 2 , K 3 加权平均得出 平均斜率K 的近似值，计算格式具有
*
yn1 yn h(1 )K1 K2 K3
K1 , K2为二阶Runge-Kutta格式的表达式 如何预报K3？
在区间[xn ,xn+q ]已知两个斜率值K1,K2，对K1,K2加权 平均得出此区间的平均斜率,从而得到y(xn+q )的预报 值yn+q
x1 x0 h; h h K1 f ( x0 , y0 ), K 2 f ( x0 , y0 K1 ) 2 2 h h K 3 f ( x0 , y0 K 2 ), f ( x0 h, y0 hK 3 ) 2 2 h y1 y 0 ( K1 2 K 2 2 K 3 K 4 ) 6
1.2649
1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.7321
改进的Euler方法
第三节
Runge-Kutta方法
拉格朗日中值定理
如果函数f ( x )满足 (1)在闭区间[a , b]上连续， )在开区间(a , b)上可导；那么, (2 至少存在一点 (a , b), 使得f (b) f (a ) f ( )(b a )成立
改进的Euler方法
综合两种方法，先用Euler法得到一个初步的近似值
yn1称为预报值，预报值yn1的精度不高,利用预报 值替代右端的yn1再直接计算，得到校正值yn1
建立预报-校正系统 预报：yn1 yn hf xn , yn h f xn , yn f xn1 , yn1 校正：yn1 yn 2
将改进Euler法推广 y n+1 =y n +h 1 K1 2 K 2 K1 f ( x n , yn ) K f ( x , y phK ) n p n 1 2
（1）
d f ( x, y) 首先希望能确定系数 1、2、p，使得到的算法格式有2阶 dx 精度，即在 yn y( xn ) 的前提假设下，使得 dy f3x ( x , y ) f y ( x , y ) Ri y( xn1 ) yn1 O(h ) dx f x ( x, y) f y ( x, y) f ( x, y) y( x )
( xn , xn1 ), st.y( xn1 ) y( xn ) y( )( xn1 xn )
y f ( x , y )
准确成立
y( xn1 ) y( xn ) hf ( , y( ))
令K f ( , y( )), 称为区间[ xn , xn1 ]上的平均斜率
第二节 改进的Euler方法
已知初值问题的一般形 式为: dy f ( x, y) dx y ( x 0 ) y0 a xb (1)
将方程两端从xn到xn+1求积分，得 y( xn1 ) y( xn )
xn1
xn
f ( x , y( x ))dx
预测—校正系统称作改进的欧拉公式。
单步显式格式
改进的Euler方法
改进Euler格式的嵌套形式
h yn1 yn f ( xn , yn )+f ( xn1 , yn +hf ( xn , yn ) 2
平均化形式： y p yn hf ( xn , yn ) yc yn hf ( xn 1 , y p ) 1 yn 1 ( y p yc ) 2 y ( x 0 ) y0
Step 1: 将 K2 在 ( xn , yn ) 点作 Taylor 展开
K 2 f ( xn ph , yn phK1 ) f ( xn , yn ) phf x ( xn , yn ) phK1 f y ( xn , yn ) O( h2 )
y( xn ) phy( xn ) O(h2 )
四、四阶Runge-Kutta方法
继续上述过程，可以进一步导出四阶Runge-Kutta 格式
yn1 yn h K1 2 K 2 2 K 3 K 4 , 6 K f x , y , n n 1 h K 2 f x n 1 , yn K1 , 2 2 K f x , y h K , n 1 n 2 3 2 2 K 4 f xn1 , yn hK 3 .
经典R-K公式 每一步计算需 要四个函数值
注意的问题
Runge-Kutta法的主要运算在于计算 Ki 的值，即计算 f 的值。计算量与可达到的最高精度阶数的关系：
每步须算Ki 的个数
2 3 4 5 6 7
n8
O(hn2 )
4 3 4 5 可达到的最高精度 O(h2 ) O(h ) O(h ) O(h ) O(h )
等于中点的斜率值f xn 1 , y( xn 1 ) , 本格式表示用中
2 2
点的斜率代替平均斜率K*
计算量：每步需要计算两次函数f的值
Q: 为获得更高的精度，应该如何进一步推广？
三、三阶Runge-Kutta方法
为进一步提高精度，设除xn+p外再考察一点
xnq xn qh,
p q 1
利用数值积分求积分项的方法离散
(1)左矩形法
xn1
ຫໍສະໝຸດ Baidu
f ( x , y( x ))dx hf ( xn , y( xn )) O( h2 )
xn
y( xn1 ) y( xn ) hf ( xn , y( xn ))
Euler格式 一阶方法
改进的Euler方法
(2)梯形法
xn1
O(h6 )
R-K(高阶)方法不唯一,选择不同的参数能得到 不同的R-K公式 R-K方法的推导是基于Taylor展开法，因而要求 解具有较好的光滑性，如果光滑性较差精度可 能不如改进Euler方法,最好采用低阶算法而将 步长h 取小。
四阶R-K方法实现
开始
输入 x 0, y 0, h, N
n N 1; x1 x0 , y1 y0
xn
h f ( x , y( x ))dx f ( xn , y( xn )) f ( xn1 , y( xn1 )) O( h3 ) 2
h yn1 yn [ f ( xn , yn ) f ( xn1 , yn1 )] 梯形格式 2
梯形格式是显式Euler格式与隐式Euler格式的 算术平均 Euler格式是显式算法，计算量小，但精度低 梯形格式，精度较高，但是隐式算法，需要通过 迭代过程求解，计算量大
斜率 一定取K1 K2 的平均值吗？
步长一定是一个h 吗？
用xn , xn1两个点的斜率值K1 , K 2的算术平均作为平均 斜率K *；xn 1处的斜率值K 2，由已知信息y n 通过Euler 公式预报
Runge-Kutta方法的设计思想 设法在[xn,xn+1]区间内多预报几个点的斜率值， 利用这些斜率值，将他们加权平均作为平均斜率 的近似，有可能构造出更高精度的计算格式
y( xn1 ) y( xn ) hK
*
寻求计算平均斜率的算法
考察欧拉法，以xn的斜率值
K1 f ( xn , yn )
作为平均斜率
考察改进的欧拉法，可以将其改写为：
1 1 yi 1 yn h K1 K 2 2 2 K1 f ( x n , yn ) K 2 f ( xn h, yn hK1 )
ynq yn qh(1 )K1 K2
K 3 f ( x n q , yn q )
因此K 3为:
利用Taylor展开法选择参数p, q, , , ，使此计算格式 具有三阶精度这类格式统称为三阶Runge-Kutta格式
R-K法的常用公式
常用的三阶R-K方法.
3 要求 Rn y( xn1 ) yn1 O(h ) ，则必须有：
1 1 2 1 , 2 p 2
这里有 3 个未知 数， 2 个方程。
存在无穷多个解。所有满足上式的格式统称为2阶 龙格 - 库塔格式。
1 (1) p 1, 1 2 为改进Euler 格式 2 1 (2)p , 1 0, 2 1为变形的Euler 格式 中点格式 2 yn h K1是Euler 法预报的中点xn 1 的近似值，K 2 近似 2 2
h yn 1 yn 6 K1 4 K 2 K 3 , K1 f x n , yn , K 2 f x 1 , yn h K1 , 2 n+ 2 K 3 f xn +1 , yn +h( K1 2 K 2 ) .
二、二阶Runge-Kutta方法
推广改进的Euler方法，在区间[xn , xn1 ]内任取一点 xn p xn ph,
* *
0 p 1
希望用xn，xn p两个点的斜率K1 , K 2加权平均得到平 均斜率K ，即 K =1 K1 2 K 2 , 1 ,2为待定参数
Y
x1 b
N 结束
例题3
例2
求解初值问题 步长h 0.1) ( 2x y y y y ( 0) 1 (0 x 1)
解
f ( x, y) y 2 x / y
y p yn hf ( xn , yn ) 改进Euler格式 yc yn hf ( xn1 , y p ) yn1 1 ( y p yc ) 2
Step 3: 将 yn+1 与 y( xn+1 ) 在 xn 点的泰勒展开作比较
yn1 yn (1 2 )h y( xn ) 2 ph2 y( xn ) O(h3 )
h2 y( xn1 ) y( xn ) hy( xn ) y( xn ) O( h3 ) 2
Step 2: 将 K2 代入第1式，得到
yn1 yn h 1 y( xn ) 2 [ y( xn ) phy( xn ) O( h2 )] yn (1 2 )h y( xn ) 2 ph2 y( xn ) O( h3 )
§2 Runge-Kutta Method
Euler法
xn Yn
0.1 1.1 0.2 1.1918
改进Euler法
0.0046 1.0959 0.0086 1.1841 0.0005 0.0009
准确解
1.0954 1.1832
|yn-y(xn)| Yn
|yn-y(xn)| y(xn)
0.3 1.2774
0.4 1.3582 0.5 1.4351 0.6 1.5090 0.7 1.5803 0.8 1.6498 0.9 1.7178 1.0 1.7848
改进Euler方法计算框图 开始
输入x0 , y0 , h, b
n 1
x1 x0 h y p y0 hf ( x0 , y0 ) yc y0 hf ( x1 , y p ) 1 y1 ( yc y p ) 2
输出x1 , y1
n n 1 x0 x1 , y0 y1
0.0125 1.2662
0.0166 1.3434 0.0209 1.4164 0.0257 1.4860 0.0311 1.5525 0.0373 1.6165 0.0445 1.6782 0.0527 1.7379
0.0013
0.0018 0.0022 0.0028 0.0033 0.0040 0.0049 0.0058
用三个点x n , x n p , x nq的斜率K 1 , K 2 , K 3 加权平均得出 平均斜率K 的近似值，计算格式具有
*
yn1 yn h(1 )K1 K2 K3
K1 , K2为二阶Runge-Kutta格式的表达式 如何预报K3？
在区间[xn ,xn+q ]已知两个斜率值K1,K2，对K1,K2加权 平均得出此区间的平均斜率,从而得到y(xn+q )的预报 值yn+q
x1 x0 h; h h K1 f ( x0 , y0 ), K 2 f ( x0 , y0 K1 ) 2 2 h h K 3 f ( x0 , y0 K 2 ), f ( x0 h, y0 hK 3 ) 2 2 h y1 y 0 ( K1 2 K 2 2 K 3 K 4 ) 6
1.2649
1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.7321
改进的Euler方法
第三节
Runge-Kutta方法
拉格朗日中值定理
如果函数f ( x )满足 (1)在闭区间[a , b]上连续， )在开区间(a , b)上可导；那么, (2 至少存在一点 (a , b), 使得f (b) f (a ) f ( )(b a )成立
改进的Euler方法
综合两种方法，先用Euler法得到一个初步的近似值
yn1称为预报值，预报值yn1的精度不高,利用预报 值替代右端的yn1再直接计算，得到校正值yn1
建立预报-校正系统 预报：yn1 yn hf xn , yn h f xn , yn f xn1 , yn1 校正：yn1 yn 2
将改进Euler法推广 y n+1 =y n +h 1 K1 2 K 2 K1 f ( x n , yn ) K f ( x , y phK ) n p n 1 2
（1）
d f ( x, y) 首先希望能确定系数 1、2、p，使得到的算法格式有2阶 dx 精度，即在 yn y( xn ) 的前提假设下，使得 dy f3x ( x , y ) f y ( x , y ) Ri y( xn1 ) yn1 O(h ) dx f x ( x, y) f y ( x, y) f ( x, y) y( x )
( xn , xn1 ), st.y( xn1 ) y( xn ) y( )( xn1 xn )
y f ( x , y )
准确成立
y( xn1 ) y( xn ) hf ( , y( ))
令K f ( , y( )), 称为区间[ xn , xn1 ]上的平均斜率