7-常微分方程数值解法

向后的欧拉方法递推公式为
yn ⇒ y ( xn )
yn +1 = yn + h ⋅ f ( xn +1 , yn +1 )
向后的欧拉方法(隐式方法 预报---校正法 向后的欧拉方法 隐式方法):预报 校正法 隐式方法 预报
1. 用欧拉方法预报 2. 用向后的欧拉方法校正 ① 用欧拉方法预报
y
(0) n +1
例题： 例题：
y ' = sin x + y y ( x0 ) = 1, x0 = 0
2.1.2 折线法（几何意义） 折线法（几何意义） 作曲线Q(t)的切 从（t0,Q0）作曲线 的切 得切线方程： 线，得切线方程：
Q
Q1
Q(t)
dQ Q = Q0 + t0 ( t − t0 ) dt = Q0 + f (Q0 , t0 ) ( t − t0 )
Qn ⇒ Q(tn )
Qn +1 = Qn + f [Qn , tn ] ⋅ h
2.1.4 欧拉法误差 利用泰勒级数得： 利用泰勒级数得：
Q ( tn +1 ) = Q(tn + h) 1 2 ≈ Q(tn ) + hQ′(tn ) + h Q′′(tn ) +K +K 2 1 2 = Q(tn ) + hf [Q(tn ), tn ] + h Q′′(tn ) + K 2
欧拉数值算法递推公式构造 2.1.1 差分法 差分法就是用差商近似代替微商， 差分法就是用差商近似代替微商，即
∆Q dQ ⇒ ∆t dt
代入微分方程得到： 代入微分方程得到：
Q (t + ∆t ) − Q (t ) = f (Q , t ) ∆t Q (t + ∆t ) = Q (t ) + f (Q, t ) ∆t
dy = f ( x, y ) x0 ≤ x ≤ xn dx y ( x0 ) = y0
数值解法就是求y(x)在某些分立的节点 xn 上的近似值 在某些分立的节点 数值解法就是求 yn，用以近似 n) 用以近似y(x
yn ≈ y ( xn )
源自文库
2、欧拉近似方法 、
2.1 简单欧拉 简单欧拉(L.Euler, 1707-1783)方法。 方法。 方法
隐式方法: 预报---校正法 隐式方法 预报 校正法
1. 用欧拉方法预报 2. 用改进的欧拉方法校正
计算公式为： 计算公式为：
y = yn + h ⋅ f ( xn , yn ) h (0) yn +1 = yn + ⋅ f ( xn , yn ) + f ( xn +1 , yn +1 ) 2
对于等间隔节点
∆t = tn +1 − tn = h tn +1 = tn + h
可以得到： 可以得到： tn Q精确值 精确值 Q近似值 近似值 t0 t1 t2
n=0,1,2K
…. …. ….
tn
…. …. ….
Q(t0) Q(t1) Q(t2) Q0 Q1 Q2
Q(tn) Qn
节点上， 在tn节点上，微分方程可以写为
= yn + h ⋅ f ( xn , yn )
用向后的欧拉方法校正---迭代 ② 用向后的欧拉方法校正 迭代
y
( k +1) n +1
= yn + h ⋅ f ( xn +1 , y )
(k ) n +1
2.2.2 改进的欧拉方法 积分法求欧拉公式时，积分采用梯形近似， 积分法求欧拉公式时，积分采用梯形近似，即得到改 进的欧拉方法。 进的欧拉方法。
2.2 改进的欧拉近似方法 对于一般的常微分方程： 对于一般的常微分方程：
dy = f ( x, y ) dx y ( x0 ) = y0
节点： x0 , x1 ,K , xi ,K , xn ,K 节点： 步长： 步长：
h = xi +1 − xi
积分法求欧拉公式时，矩形法采用前一点函数值求积， 积分法求欧拉公式时，矩形法采用前一点函数值求积， 若利用后一点，即为向后的欧拉方法。 若利用后一点，即为向后的欧拉方法。
t1
= Q0 + f [Q(t0 ), t0 ] ⋅ h
t0
同样， 积分采用矩形近似， 同样，在[t0,t2] ，积分采用矩形近似，得：
Q(t2 ) = Q0 + ∫ f [Q(t ), t ] ⋅ dt
t2 t0
= Q0 + ∫ f [Q(t ), t ] ⋅ dt + ∫ f [Q(t ), t ] ⋅ dt dt
即
∫
y y0
dy = ∫ f ( x, y )dx
x0 xn+1 xn
x
y ( x) = y0 + ∫ f ( x, y )dx
x0
x
y ( xn +1 ) = y ( xn ) + ∫
f ( x, y )dx
= y ( xn ) + h ⋅ f xn +1 , y ( xn +1 )
(0) n +1
例题： 例题： y ' = sin x + y
y ( x0 ) = 1, x0 = 0
3、龙格-库塔 、龙格 库塔 库塔(R-K)方法 方法
几种方法的比较 对于一般的常微分方程： 对于一般的常微分方程：
dy = f ( x, y ) dx y ( x0 ) = y0
y ( xn +1 ) = y ( xn ) + ∫
xn+1
xn
f ( x, y )dx
h = y ( xn ) + ⋅ f xn , y ( xn ) + f xn +1 , y ( xn +1 ) 2
{
}
改进的欧拉方法递推公式为
h yn +1 = yn + ⋅ [ f ( xn , yn ) + f ( xn +1 , yn +1 )] 2
∫
Q( t )
Q0
dQ = ∫ f Q ( t ) , t dt t0
t t
Q ( t ) = Q0 + ∫ f Q ( t ) , t dt t0
积分采用矩形近似， 在[t0,t1] ，积分采用矩形近似，得：
Q(t1 ) = Q0 + ∫ f [Q(t ), t ] ⋅ dt
dQ = f (Q, t ) dt Q (t0 ) = Q0
欧拉数值算法就是由初值通过递推求解，递推求解 就是从初值开始，后一个函数值由前一个函数值得到。 就是从初值开始，后一个函数值由前一个函数值得到。关 键是构造递推公式。 键是构造递推公式。
Q0 ⇒ Q1 ⇒ Q2 ⇒ K ⇒ Qn
Q1
Q(t)
Q2
dQ Q = Q1 + t1 ( t − t1 ) dt Q1 + f (Q1 , t1 ) ( t − t1 )
交点纵坐标为： 与t=t2交点纵坐标为：
t0 t1 t2
t3 t4 t5 t6
t
Q1 ⇒ Q(t1 )
Q2 = Q1 + f (Q1 , t1 ) ⋅ h
按照相似的方法， 按照相似的方法，从（tn,Qn） Q 作曲线Q(t)在(tn,Q(tn))的切 作曲线 在 的切 线的平行线，得直线方程： 线的平行线，得直线方程：
Q1
Q(t)
Q2
Q
Qn + f (Qn , tn ) ( t − tn )
交点纵坐标为： 与t=tn+1交点纵坐标为：
t0 t1 t2 t3 t4 t5 t6 t
Qn +1 = Qn + f (Qn , tn ) ⋅ h
折线近似曲线， 增大 增大， 折线近似曲线，n增大，误差变大
2.1.3 数值积分 积分得： 对微分方程作从 t0 到 t 积分得：
科学计算与MATLAB 科学计算与MATLAB
中南大学材料科学与工程学院 2011.10
第七讲
常微分方程数值解法
内容提要
引言 欧拉近似方法 龙格-库塔 龙格 库塔(R-K)方法 方法 库塔 MATLAB的常微分方程函数 的常微分方程函数 小结
1、引言 、
物理学所研究的各种物质运动中， 物理学所研究的各种物质运动中，有许多物质运动的 过程是用常微分方程来描述的。 过程是用常微分方程来描述的。 例如,质点的加速运动，简谐振动等。 例如,质点的加速运动，简谐振动等。
Q3 = Q2 + f (Q2 , t2 ) ⋅ h M M M M M M
简单欧拉方法程序
function [outx,outy]=MyEuler(fun,x0,xt,y0,PointNum) %MyEuler 用前向差分的欧拉方法解微分方程 %fun 表示 （x，y） 表示f（ ， ） %x0,xt表示自变量的初值和终值 表示自变量的初值和终值 %y0表示函数在 处的值 其可以为向量形式 表示函数在x0处的值 表示函数在 处的值,其可以为向量形式 %PointNum表示自变量在 表示自变量在[x0,xt]上取的点数 表示自变量在 上取的点数 if nargin<5 | PointNum<=0 %如果函数仅输入 个参数值，则PointNum默认值为 如果函数仅输入4个参数值 默认值为100 如果函数仅输入 个参数值， 默认值为 PointNum=100; end if nargin<4 %y0默认值为 默认值为0 默认值为 y0=0; end h=(xt-x0)/PointNum;%计算步长 计算步长h 计算步长 x=x0+[0:PointNum]'*h;%自变量数组 自变量数组 y(1,:) = y0(:)';%将输入存为行向量，输入为列向量形式 将输入存为行向量， 将输入存为行向量 for k = 1:PointNum f=feval(fun,x(k),y(k,:));%计算 计算f(x,y)在每个迭代点的值 计算 在每个迭代点的值 f=f(:)'; y(k + 1,:) =y(k,:) +h*f; %对于所取的点 迭代计算 值 对于所取的点x迭代计算 对于所取的点 迭代计算y值 end outy=y; outx=x; %plot(x,y)%画出方程解的函数图 画出方程解的函数图
Q(tn +1 ) = Q(tn ) + f [Q(tn ) , tn ] ⋅ h
作如下近似： 作如下近似：
Qn ⇒ Q(tn )
则得到欧拉解法递推公式的一般形式： 则得到欧拉解法递推公式的一般形式：
Qn +1 = Qn + f (Qn , tn ) ⋅ h
具体求解过程为： 具体求解过程为：
Q1 = Q0 + f (Q0 , t0 ) ⋅ h Q2 = Q1 + f (Q1 , t1 ) ⋅ h
作如下近似
Qn ⇒ Q(tn )
Qn +1 = Qn + f [Qn , tn ] ⋅ h
局部截断误差 由泰勒展式可以看出，在欧拉方法中，假定Q 由泰勒展式可以看出，在欧拉方法中，假定 n 是精 确的， 计算Q 所产生的误差的数量级为O(h2) 。 确的，由Qn计算 n+1所产生的误差的数量级为 整体截断误差 如果再考虑到Q 本身的误差， 如果再考虑到 n本身的误差， 计及误差的累积效果 和局部截断误差，这种误差称为欧拉方法整体截断误差， 和局部截断误差，这种误差称为欧拉方法整体截断误差， 数量级为O(h) 。 数量级为 欧拉法是数值方法解微分方程最简单的一种， 欧拉法是数值方法解微分方程最简单的一种，精度 不高，实际应用不多， 不高，实际应用不多，但它可以反映数值求解微分方程 的特征。 的特征。
dv F =m dt
d x 2 2 +ω x = 0 2 dt
简单问题可以求得解析解，多数实际问题靠数值求解。 简单问题可以求得解析解，多数实际问题靠数值求解。
2
一阶常微分方程(ODE )初值问题 ： 一阶常微分方程 初值问题 ODE ：Ordinary Differential Equation
t1 t2
= Q(t1 ) + f [Q(t1 ), t1 ] ⋅ h
t0
t1
同样， 积分采用矩形近似， 同样，在[t0,tn+1] ，积分采用矩形近似，得：
Q(tn +1 ) = Q0 + ∫
tn+1
= Q(tn ) + f [Q(tn ), tn ] ⋅ h
作如下近似 得：
t0
f [Q(t ), t ] ⋅ dt
t0 t1 t2
t3 t4 t5 t6
t
Q1 = Q0 + f (Q0 , t0 ) ( t1 − t0 ) =Q0 + f (Q0 , t0 ) ⋅ h
此切线与t=t 交点纵坐标为： 此切线与 1交点纵坐标为：
Q0 ⇒ Q (t0 )
作曲线Q(t)在 Q 从（t1,Q1）作曲线 在 (t1,Q(t1))的切线的平行线， 的切线的平行线， 的切线的平行线 得直线方程： 得直线方程：