第四章常微分方程的数值解法
Numerical Solutions to Ordinary
Differential Equations
概述
(())
(())
Ι,Problem
00
dy
=f x,y
dx
yx =y
12 1 2 1 2
L,y,y,|f(x,s,y)-f(x,y)| L|y |t,-y ≤
!一阶常微分方程初值问题:
! f(x,y)在[a,b]上连续,且满足Lipschitz条件:
!实际工程技术 生产 科研上会出现大量的微分方程问题很难得到其解析解 有的甚至无法用解析表达式来表示因此只能依赖于数值方法去获得微分方程的数值解求函数y=y(x)满足
y?(x) =f(x,y(x))
则初值问题Problem I有唯一解
!微分方程的数值解法
"不求y=y(x)的精确表达式,而求离散点x
0
,x
1
,…x
n
处的函数值
"设Problem I的解y(x)的存在区间是[a,b] 初始点x
0
=a 取
[a,b]内的一系列节点x
0
,x
1
,...,x
n
a= x
0
< x
1
<…< x
n
=b
一般采用等距
a b
x
0
x
1
x
2
..,x
n-1
x
n
"用数值方法 求得y(x)在 节点x
k
的值y(x
k
) 的
值 用y
k
表示 y
k
y(x
k
) y
0
,y
1
,...,y
n
微分方程的数值解
"求y(x)————>求y
0
,y
1
,...,y
n
! 程
00
n+1 n,n n-1 n-2 n-m
y=y(x)
y=g(h,x y,y,y,...y )
!!方法 采用 式 法
[a,b]n等分 a= x
0
< x
1
<…< x
n
=b h= x
k
=a+kh
b-a
n
10 112 nm nm n n
yy yy yyy
++
→ → →→ → →→ → →
! 式
!Taylor 式
!数值 分法
!,用 微,
≈
n+1 n
nn+1n
y(x )-y(x )
y'(x ) h=x -x
h
00
n+1 n n n
y(x )=y
n=0,1,2,
y=y+hf(x,y)
——? 式 Euler Scheme
4.1? 式
∴≈
n+1 n
n n
y(x )-
y(x
y(x )
f(,)x),
h
(())
(())
'
00
y=f x,y
yx =y
∴≈
n+ n n1 n
y(x )y( +hf(x,) y(xx ))
∴
nnn+1 n
=y +hf(x,yy )
￠￡?￥
1,y(x) 点P
0
(x
0
,y
0
)且在点(x,y)的§currency1
',f(x,y),
2,y(x)在点P
0
(x
0
,y
0
)的§
currency1方程
3,? 2 P
1
f (x
1
,y
1
) '“?'“?y(x)的§currency1 在其上取点
P
2
(x
2
,y
2
) 依此? …
4.ficurrency1P
0
P
1
P
2
P
n
flcurrency1y(x)的
——?ficurrency1法
x
p
0
p
1
p
2
p
3
p
4
x
0
x
1
x
2
x
3
x
4
y
y(x)
(())
(())
'
00
y=f x,y
yx =y
000 0
y=y +f(x,y )(x-x )
在§currency1上取点P
1
(x
1
,y
1
)
01 00
=y + f(xy,yh )
y
1
是Euler 式–求
!用? 微
≈
n+1 n
n+1 n n+1 n+1 n
y(x )-y(x )
y'(x ) y[x,x ]= h=x -x
h
00
n+1 n n+1 n+1
y(x )=y
n=0,1,2,
y=y+hf(x,y)
——?式? 式
法 续
∴≈
n+
n+
1
1n
n+1
y(x )-y(x )
f(x,)y(x )
h
(())
(())
'
00
y=f x,y
yx =y
≈?
n+1 nn1 +n+ 1
+hf(x,y(x ) y() xy(x ))
∴
nnn1 +1+ n+1
=y +hf( yx,y )
here
· 用?式?法 一?
解方程解出y
n+1
的?式表达式,其?…‰
!用? 微
≈
n-1 n+
n+1 n-1
1n
y(x )-y(x )
y'(x ) y[,]=xx
2h
≈
n+1 n-1
n n
y(x )-y(x )
f(x y(,)x)
2h
n-1
00
n+1 nn
y(x )=y
n=0,1,2,
=+2hf(x )y,yy
—— `?法
· ′?用
法求出y
1
用
`?法
法 续
n+1 n n n
y=y+hf(x,y)
式? ˉ?式ˉ
式
n+1 n n+1 n+1
y=y+hf(x,y)
式
n+1 n-1 n n
y=y+2hf(x,y) `?式
˙¨ y(x
n+1
)-y
n+1
˙¨
…￥1?设yn=y(xn),第n 是精确的?,
Rn+1=y(xn+1)-yn+1?法的˙¨,
…法的精ˇ是一阶
…￥2 — 法的˙¨ O(h
p+1
),法有p阶精ˇ.
·
无yn=y(xn)? Rn+1?˙¨
分析 其˙¨ O(h
2
),Taylor 式式分析分析
Euler 式
nn+1 nn
=y +hf(x,yy )
yn=y(xn) y(xn+1)-yn+1 =O(h
2
)
∵
n+1 nnn
=+hf(xy,yy )
nnn
=+hf(y(x ) y,(xx ))
y?(x) =f(x,y(x))
'
=+hy(x) y(x)
T
'
aylor
'
n+1 n n n n+1
2
n
'
y(x )=y(x +h) === y(x )+hy (x )+,(x,x
y)
2
)
(
h ξ
ξ
∈
!
∵
''
2
n+1
() 2
2
n+1
=y( hx)-y h)=O(
y ξ
∴
`?法的˙¨
…￥3?设yn=y(xn),y
n 1
=y(x
n 1
),Rn+1=y(xn+1)-yn+1
`?法的˙¨,
n-1 nn n+1
=+2y hf(x,yy)
…式?法的精ˇ是一阶 `?法的精ˇ是`阶
`?法 其˙¨
yn=y(xn),y
n 1
=y(x
n 1
),
nnn-1
=+2hfy(x ) yx,(x( ))
'
n-1 n
=+y( 2hx) y(x)
Taylor
'
2
''
n+1 n n n n n+
'
1
''
3
n
y(x )=y(x +h)=== y(x )+hy (x )
h
y(x )
2!
++,(
y(
)
)
x,xh
ξ
ξ ∈
3!
∵
n+1 n+
3
1
y(x )-y =O(h )∴
n-1 nn n+1
=+2y hf(x,yy)
'
n-1 n n n
2
n-1
''
3
n n
'
''
(-h)
y(x )
2
y(
y(x )=y(
)
x -h)=y(x )-hy (x )+ +,(x,x )(-h)
!
η
η ∈
3!
上?式 a? ′
''' '''
'
n-1
() ()
3
3
+ nn1
-y(x ) 2hy (x=y(x ) )+ h
yyξη+
!
∴
`?法的˙¨ O(h
3
) 其精ˇ是`阶
,求 的 值
2
(0) 1
dy x
y
dx y
y
=?
=
0.1,0.2,1.0,x = "
解,?
2
(,)
x
fxy y
y
=?
00
0,1,0.1xyh===
式 得:
1
(,)
nn nn
yyhf x y
+
=+
0
1y =
10 00
0
(,)10.1(1 )1.1
1
yyhfxy=+ =+×?=
21 11
20.1
(,) 1.1 0.1 (1.1 ) 1.191818
1.1
yyhfxy
×
=+ =+×? =
32 22
(,) 1.277438yyhfxy=+ =
……
其精确解 21yx=+
¨?o 其?分析
111
()
kkk
e y x y
+++
=?
表? 出¨ 在?
10 9 9 9
( ) (,( )) 1.7330815yyxfxyx=+ =
#
˙¨
10 10
( ) 0.00103yx y?=
#
而¨ 是 10 10
( ) 0.05272yx y?=
数值 分法
(())
(())
'
00
yx
yy
=y
=f x,
n+1 n+1
nn
x
'
x
xx
y( ) f(x,y)xdx dx∴=
∫∫
n+1
n
x
x
1
y(f(x,)(( x)))
nn
yx yx dx
+
=?
∫
a 的… 分用数值 分 式求 值
1 用 数值 分 式
n+1
n
1
x
x
(f(x,) (,)y(x) y() )
nnn n
dfxx xxx
+
≈
∫
1
y( )-y( ) ( y(,))
nn nn
hf xxx x
+
≈∴
1
y( )y( ) + (,)y( )
nnnn
xhfxx
+
≈
+(,)
1
yyy
n
hf x
n nn
↓↓ ↓
=
+
2 用 式
1
1
n+
n
nn
x
nn+11
(
+
)
2
x
[y(xf(x,) f(x,))yy(x)+f(x,)]
nn
xx
dx
+
≈
∫
n+1 nn+
2
1nn+1n
yf(x,) fy(x )-y (x,(x ) y(x(x ) [ )])+
h
∴≈
nnn +1n+ +
2
1 1
y=y [f(x +],) f(x,)yy
h
∴+
—— 式
n+1 n n+1 n+1
y=y+hf(x,y)
n+1 n n n
y=y+hf(x,y)
# 式,?示? 式,?式? 式? 得
# 式是?式? 式,其精ˇ `阶
式的精ˇ
y
n
=y(x
n
),Talor 式有分析 式是?式 式 其˙¨ O(h
3
)
用 用到到`?函数的Taylor 式
…? 式 的精ˇ是2阶的.
n nn+1n1
2
+1+
yf=y [ (x,) f(x,)y+]y
h
+
f(x
n+1
,y
n+1
)=f(x
n+1
,y(x
n+1
)+(y
n+1
-y(x
n+1
))
=f(x
n+1
,y(x
n+1
))+f
y
(x
n+1
,)(y
n+1
-y(x
n+1
)) (x
n
x
n+1
)
=y’(x
n+1
)+f
y
(x
n+1
,)(y
n+1
-y(x
n+1
))
=y’(x
n
)+hy”(x
n
)+O(h
2
) +f
y
(x
n+1
,)(y
n+1
-y(x
n+1
))
=f(x
n
,y
n
)+hy”(x
n
) +f
y
(x
n+1
,)(y
n+1
-y(x
n+1
)) +O(h
2
)
y(x
n+1
)= y(x
n
+h) = y(x
n
)+hy’(x
n
)+h
2
y”(x
n
) /2 +O(h
3
)
=y
n
+hf(x
n
,y
n
)+h
2
y”(x
n
)/2+O(h
3
)
=y
n
+hf(x
n
,y
n
)/2+h[f(x
n
,y
n
)+hy”(x
n
)]/2+O(h
3
)
而y(x
n+1
)=y
n+1
+h f
y
(x
n+1
,)[y(x
n+1
)-y
n+1
] /2 +O(h
3
)
y(x
n+1
)-y
n+1
=h f
y
(x
n+1
,)[y(x
n+1
)-y
n+1
] /2 +O(h
3
)
y(x
n+1
)-y
n+1
=O(h
3
)/[1-hf
y
(x
n+1
,)/2]=O(h
3
)
式的?˙¨ O(h
3
) 其精ˇ是精ˇ是2阶
f(x
n+1
,y
n+1
)=f(x
n
,y
n
)+hy”(x
n
) +f
y
(x
n+1
,)(y
n+1
-y(x
n+1
)) +O(h
2
)
y(x
n+1
)= y
n
+hf(x
n
,y
n
)/2+h[f(x
n
,y
n
)+hy”(x
n
)]/2+O(h
3
)
= y
n
+hf(x
n
,y
n
)/2+h [f(x
n+1
,y
n+1
)
-f
y
(x
n+1
,)(y
n+1
-y(x
n+1
)) +O(h
2
)]/2+O(h
3
)
= y
n
+h[f(x
n
,y
n
) + f(x
n+1
,y
n+1
) ]/2
+ h f
y
(x
n+1
,)(y(x
n+1
) - y
n+1
) /2 +O(h
3
)
式的?用
dy
=y,y(0)=1.
dx
4.1 用 式求初值问题的解在x=0.01上的值y(0.01).
10 01
h
y=y+ [y+y]
2
解,取h=0.01 x
0
=0 y
0
=y(0)=1 则 y(0.01) y
1
f(x,y)=y 式,
≈+ + =
23
0.01
0.01 0.01
y(0.01)=e =1+0.01+ + +...
2! 3!
.
..
!
2
001
1 0 01 1 01005
2
n+1 n n n n+1 n+1 n n n+1
hh
y =y + [f(x,y )+f(x,y )]=y + [y +y ]
22
h
2
h
0
2
1
1
=
-
1+
yy
≈
10
2
2
2
h
y=(1+ )( )y
2
hh
(1+ ) + =1.01005
2
hh
1+ + +...
24
4
于?…‰
解析解
x
y=e
式的?
精ˇ有–,”?式,? 求解,
量大
式法
式,” 初值,且第2 初值只能 其
方法给出,能?o的 精ˇ有影响
`?法
…‰?式?法最简?,精ˇ低
法
4.2 改 的的Euler法
!Euler 式 y
n+1
=y
n
+hf(x
n,
y
n
)?式 一阶
! 式 y
n+1
=y
n
+h[f(x
n
,y
n
)+f(x
n+1
,y
n+1
)]/2?式 `阶
!Euler 式 量小 精ˇ低
! 式 量大 ˇ
综合? 式 出预报—校 式
n
_
n+1 nn
=y +hf(x,yy )
预报
h
nnnn
_
n
n+12
+1
+1
=y + [f(x,y )+f(x,]y )
y
校
——改 的Euler法写成嵌套 式 n+1
y
h
nnnn+12 nnn
=y + [f(x,y )+f(x y+hf(x,y,))]
式? 法
化 式
n+1
y
nnn nnn+1 n
1
2 n
y+y+hf(x,y) y+hf(=[ x,yhf(x,+ ))]
pn nn
cn n+1p
1
n+1 p c
2
y=y+hf(x+y)
y=y+hf(x +y)
y=(y+y)
4.4 用改 的Euler法解初值问题在区间[0,0.4]上
h=0.1的解 并? 与精确解的 异
2
dx
=y
dy
y(0) = 1
说 精确解y=1/(1-x)
解 Euler法的具 式 y
n+1
=y
n
+hy
n
2
改 的Euler法的具 式
2
pn n
2
cn p
1
n+1 p c
2
y=y+hy
y=y+hy
y=(y+y)
x
0
=0,h=0.1,则
x
1
=0.1,x
2
=0.2,x
3
=0.3,x
4
=0.4
y
1:
y
p
=y
0
+0.1y
0
2
=1+0.1·1
2
=1.1
y
c
=1+0.1X1.1
2
=1.121
y
1
=(1.1+1.121)/2 1.1118
求y
2,
y
3,
y
4
见P95表:
y(x
n
)-y
n
随n 大而 大 表 ¨?
·
1 y(x
n
)=y
n
导改 的Euler法的˙¨
O(h
3
),具有`阶精ˇ
2 改 的Euler法也 写成如 化 式
n+1 n 1 2
1nn
2n+1,n1
1
y=y+(k+k)
2
k=hf(x,y)
k=hf(x y+k)
n+1
y
h
nnnn+12 nnn
=y + [f(x,y )+f(x y+hf(x,y,))]
HW,
p.117 #3 4
龙格—库塔方法
1n+1 n
1nn
y=y+
k=hf(x,y)
k
改 的Euler 式,
1
12n+1 n
1nn
2n+1
2
,n 1
y=y+
k=hf(x,y)
k=hf(x y+k)
(k +k )
'(,)
()
00
yfxy
yx y
=
=
Euler 式 y
n+1
=y
n
+hf(x
n,
y
n
)
写成精ˇ 一阶精ˇ ``阶阶
n+1 n
(x ) (x )
y-y
h
Lagrange?值…?
ξ?∈
nn+1
(x,x )
ξ
'
=y( )
ξ∴
n+1 n
'
y(x ) = y( + hy) (x )
ξξξ而
'
=f((),y(y ))
n+1
*
n
y(x ) = (x ) +kyξξ∴
n+1 n
hf(y(x ) = y,y()+ (x ))
y(x)在[x
n
,x
n+1
]上的? '“
*
nn 1
k=hf(x,y)=k
!取
—Euler 式 !取
12
kk
*
22
k= +
—改 Euler 式
Euler 式用一 点的值k
1
k
*
的 值 而改 的
Euler 式用` 点的值k
1
k
2
的? 值? k
*
值 改
的Euler法?Euler法精ˇ
龙格 库塔法
Runge-Kutta法的 在 x
n,
x
n+1
内多预报￠ 点的k
i
值并用其 权? 值? k
*
值 而构造出具有更 精ˇ的
式
4.3.2 `阶龙格—库塔方法
k
1
与k
2
的 权? 来 k
*
设
αβ
n+1 n 1 2
1n
12
n
2n n1
y=y+k+k
k = hf(x,y ) (*)
k=hf(x+h,y+ k)
ww
其?w
1,
w
2,
α,β 待…参数 适当选取参数 使(*)式的精ˇ `阶 使其˙¨ O(h
3
)
y(x
n
)=y
n
泰勒 式
2
''' 3
h
n+1n n n2
y(x )=y(x )+hy(x )+ y (x )+O(h )
∵
'
nnn nn
y (x ) = f(x,y(x )) = f(x,y )
'' ' ' ' '
xxy xy
( y (x)=(y (x)) =[f(x,y)] =f +f y (x)=f +f f(x,y) )∵
''
nxnnynn nn
y (x )=f (x,y )+f (x,y )f(x,y )
∴
2
3
h
nnn xnn+1 n ynn nn2
=y+hf(x,y)+ [f(x,y)y(x +f (x,y )f(x,y )]+O(h) )
`阶龙格 库塔法
αβα β
2
2 n 1 nn xnn 1ynn
k =hf(x + h,yn+ k )=h{f(x,y )+ hf (x,y )+ k f (x,y )+O(h )}
多?函数的泰勒 式
αβ
n+1 n 1 2
1n
12
n
2n n1
y=y+k+k
k = hf(x,y ) (*)
k=hf(x+h,y+ k)
ww
α
β
∴
2
n1nn 2nn 2xnn
23
2ynnn
+1
n
n
=y +w hf(x,y )+w hf(x,y )+w h f (x,y )
+w h f (x,y )f(x,y )+O(h
y
)
α
β
2
h
nnxn2n+ 2
3
ynn nn2
1 12
=y+ hf(x,y)+ [ f(x,y)
+ f (x,y )f(x,y )]+O(h )
2
(w +w )
w
2w
y
(2)
2
3
h
nnn xnnynnnn2n+1
=y +hf(x,y )+ [f (x,y )+f (x,y )f(x,y )]+O(h ) y(x ) (1)
则有
w
1
+w
2
=1
w
2
α=1/2
w
2
=1/2
· 上述方程组有四 未知量 只有三 方程 有无穷多组解
<1>与<2> 使 y(x
n+1
)-y
n+1
=O(h
3
)
`阶龙格 库塔法
!取一组解便得一种`阶龙 库 式
αβ
n+1 n 1 2
1nn
2
12
nn1
y=y+k+k
k=hf(x,y)
k=hf(x
w
+h,y+k)
w
w
1
+w
2
=1
w
2
α=1/2
w
2
=1/2
!当w
1
=w
2
=1/2,= =1′`阶Runge-Kutta 式
y
n+1
=y
n
+k
1
/2+k
2
/2
k
1
=hf(x
n
,y
n
)
k
2
=hf(x
n
+h,y
n
+k
1
)
此 改 得Euler法
!取w
2
=0 w
2
=1 =1/2 =1/2
y
n+1
=y
n
+k
2
k
1
=hf(x
n
,y
n
)
k
2
=hf(x
n
+h/2,y
n
+k
1
/2)
此?点法?变 的 Euler 式三阶龙格 库塔法
!三阶龙格 库塔法是用三 值k
1
,k
2
,k
3
的 权? 来 k*
有
y
n+1
=y
n
+c
1
k
1
+c
2
k
2
+c
3
k
3
k
1
=hf(x
n
,y
n
)
k
2
=hf(x
n
+a
2
h,y
n
+b
21
k
1
)
k
3
=hf(x
n
+a
3
h,y
n
+b
31
k
1
+b
32
k
2
)
! 使其具有三阶精ˇ 必须使˙¨ O(h
4
)
!? `阶龙格 库塔法的 导 c
1
,c
2
,c
3
,a
2
,a
3
,b
21
,b
31
,b
32
满足
c
1
+c
2
+c
3
=1
a
2
=b
21
a
3
=b
31
+b
32
c
2
a
2
+c
3
a
3
=1/2
c
2
a
2
2
+c
3
a
3
2
=1/3
c
3
a
3
2
=1/6
方程组解
得三阶龙格-库塔 式
Kutta 式
k
n+1
=y
n
+(k
1
+4k
2
+k
3
)/6
k
1
=hf(x
n
,y
n
)
k
2
=hf(x
n
+h/2,y
n
+k
1
/2)
k
3
=hf(x
n
+h,y
n
-k
1
+2k
2
)
四阶龙格 库塔法
!? 出四阶龙格-库塔 式 常用的有
经典Runge-Kutta法
y
n+1
=y
n
+(k
1
+2k
2
+2k
3
+k
4
)/6
k
1
=hf(x
n
,y
n
)
k
2
=hf(x
n
+h/2,y
n
+k
1
/2)
k
3
=hf(x
n
+h/2,y
n
+k
2
/2)
k
4
=hf(x
n
+h,y
n
+k
3
)
˙¨ O(h
5
)
!还有,Gill 式及m (m>4)阶龙格 库塔法
m>4′ 量太大 精确ˇ不一… 有′会降低龙格 库塔法
于经典的四阶Runge-Kutta法给出如? 法求解 dy/dx=f(x,y) a x b
y (a)=y
0
!Step 1,输入a,b,y
0
及N
!Step 2,(b-a)/N=>h,a=>x,y
0
=>y
!Step 3,输出 (x,y)
!Step 4,For I=1 T0 N
" hf(x,y)=>k
1
" hf(x+h/2,y+ k
1
/2)=> k
2
" hf(x+h/2,y+k
2
/2)=>k
3
" hf(x+h,y k
3
)=>k
4
" y+(k
1
+2k
2
+2k
3
+k
4
)/6=>y
" x+h=>x
" 输出(x,y)
!END
龙格 库塔法
用四阶经典Runge Kutta方法解初值问题
()
=
=
10
2
y
y
x
y
dx
dy
2.0=h
(1)求
1
y
2.0,1,0
00
=== hyx
() 1
2
,
0
0
0001
=?==
y
x
yyxfK
9181818.0
2
2
2
2
2
,
2
10
0
10
1002
=
+
+
+=
++=
K
h
y
h
x
K
h
y
K
h
y
h
xfK
9086375.0
2
2
2
2
2
,
2
20
0
20
2003
=
+
+
+=
++=
K
h
y
h
x
K
h
y
K
h
y
h
xfK
()
()
8432399.0
2
,
30
0
30
3004
=
+
+
+=
++=
hKy
hx
hKy
hKyhxfK
()
1832293.1
2
6
432101
=
++++= KKKK
h
yy
() 1832160.14.112
11
==+= xxy
5
1
103.1
×=e
龙格 库塔法
(2)求
2
y
1
0.2,0.2xh==
1832293.1
1
=y
() 8451714.0
2
,
1
1
1111
=?==
y
x
yyxfK
()
2111
22
,0.7946656
hh
Kfx yK=++=
()
3112
22
,0.7874945
hh
Kfx yK=++=
()
4113
,0.7440375KfxhyhK=++=
()
21 1 2 3 4
2 1.3416803
6
h
yy KK KK=+ ++ + =
() 3416408.112
22
=+= xxy
5
2
100.4
×=e
变 龙格 库塔法
!用户 出问题I
问题 如￡判˙|y(x
n
)-y
n
|<ε
精确值y(x
n
)未知如￡取h=?
解 如用p阶龙格 库塔法˙¨ O(h
p+1
)
y’=f(x,y)
y(x
0
)=y
0
求¨ < 10
-8
求数值解
x
n
h/2 x
n+1
h
!如 x
n
-----------------$x
n+1
y
n
=y(x
n
)y
n+1
(h)
则 y(x
n+1
)-y
n+1
(h)
ch
p+1
fi半x
n
$x
n+h/2
$x
n+1
分? y(x
n+1
)的 值y
n+1
(h/2)
则 y(x
n+1
)-y
n+1
(h/2)
2c(h/2)
p+1
h
2
n+1
()
n+1
(h) p
n+1 n+1
y(x )-y
1
y(x )-y 2
∴≈
h
2
h
2
nn +1+1
( ()
(h)
)
n+1
p
n +1
1
[]
2
y-yyy
1
(x )- ≈
$
变 龙格 库塔法
…? 于问题I—用P阶龙格 库塔法 y(x
n+1
)在 fi半
的 值分别 y
n+1
(h)
,y
n+1
(h/2)
则有¨ 式
p
h
2
n+1
h
2
n+1
1
2
()
n
()
(h)
n+
1
1 1+
y(x || )-y yy|- =|
≈?
· 1
0
¨ 的事?估 法
2
0
停机准则 <ε ( 保 |y(x
n+1
)-y
n+1
(h/2)
|< )
解
h取大˙¨ ch
p+1
大 不精确
h取小 运 量大 多 舍入¨?大解决策略 变 龙格 库塔法
If > fi半反复 直至 < 止取h 最?一次的 y
n+1
最?一次 的结果
Else if ( < ) 倍反复 直至 > 止最?一次运 的 一次 结果 –?
HW,p.117 #4
用经典四阶Runge-Kuttta法
4.5 收敛与?…‰
于常微分方程初值问题
!求y=y(x)--------$求y(x
n
)------------$求y
n
x
n
=x
0
+nh
离散化 种 式
Euler法 改 的Euler法 龙格 库塔法 是? 法数值解法
(())
(())
'
00
y=f x,y
yx =y
法 y
n+1
′只用到 一 的结果y
n
!?式? 法
y
n+1
=y
n
+hφ(x
n
,y
n
,h)
(x,y,h) 量函数 它依赖于f,仅是x
n
,y
n
,h的函数
!Def:— 数值方法 于固…的x
n
=x
0
+nh,当
h->0(n->∞)′ y
n
->y(x
n
) 则 法是收敛的
0
(0) )(
lim
n
h
n
n
yx y
→
→∞
=
(x
n
=x
0
+nh 固…值)
收敛与?…‰
· 数值方法是ˉ收敛取决于¨ y(x
n
)-y
n
的变化情况
于p阶 式,其˙¨ O(h
p+1
)
y(x
n+1
)-y
n+1
=O(h
p+1
) 其?…y
n
=y(x
n
).
虽h->0′˙¨ y(x
n+1
)-y
n+1
->0 并不能说
法收敛 因其 y
n
=y(x
n
)不成 此 入——
Def,y(x
n
)-y
n
法的解y
n
的?˙¨
˙¨?˙¨ 的区别是
法收敛%
0
(0) )(
lim
n
h
n
n
yx y
→
→∞
=
h->0′?˙¨ y(x
n+1
)-y
n+1
->0
Φ
ΦΦ≤
__ __
|(x,y,h)- (x,y,h)|L|y-y|
Th —? 法y
n+1
y
n
+h φ(x
n
,y
n
,h)具有p阶精ˇ 且 量函数
φ(x,y,h) 于y满足Lipschitz条件
—初值y
0
是准确的(y
0
=y(x
0
)) 则其?˙¨
y(x
n
)-y
n
=O(h
p
).
·
—? 法满足?上条件,? 其收敛改 Euler法的收敛‰
!判…改 Euler法的收敛‰
y
n+1
=y
n
+h[f(x
n
,y
n
)+f(x
n
+h,y
n
+hf(x
n
,y
n
))]/2
则 (x,y,h)=[f(x,y)+f(x+h,y+hf(x,y))]/2
— y
o
=y(x
0
)
≤
__ __
|f(x,y)-f(x,y)| L|y-y|
f(x,y) 于y满足 --条件:
1
2
ΦΦ
≤
__ __ _
__
_
|(x,y,h)- (x,y,h)|
[|f(x,y)- f(x,y)|+|f(x+h,)- f(x+y+hf(x,y) y+hh,f(x,y))|]
则
1
2
≤
__ __ __
[L|y- y |+L|y+hf(x,y)- y-hf(x,y)|]
≤
__ __ __ __
2
1
[L|y- y |+L|y- y |+L h
h
L(1+ L)
2
|y-y|]= |y-y|
2
…h h
0
则
ΦΦ≤
__
0
__
|(x,y,h)- (x,y,h
h
L(1+ L)
2
)| |y-y|
(x,y,h) 满足 条件,Th改 的Euler法收敛
! —f(x,y) 于y满足 条件,且 y
0
=y(x
0
)′
Euler法 准四阶龙格 库法也收敛
…‰
! 在 收敛‰′ 一般 数值方法 程是准确的 实际上 并 如此初始值y
0
有¨ δ
0
y
0
-y(x
0
).
的 一 有舍入¨
! 初始¨ 在 程的? 是 的 还是
‰ 是?…‰问题
Def 设在节点x
n
处用数值法得到的? 数值解 y
n
而实际
得到的 值 第 第n 数值解的
n
~
y
δ
nnn
~
=y -y
Def —一种数值方法在节点x
n
处的数值解y
n
的 δ
n
0
而在的 节点值y
m
(m>n)上有 δ
m
.
当|δ
m
| |δ
n
| (m=n+1,n+2,……) 则 数值 法是?…的
! 分析 法的数值?…‰ 常 方程
y’=λy (λ<0)
Euler法的?…‰
Euler法 y
n+1
=y
n
+hf(x
n
,y
n
)
方程 y’=λy (λ<0)
y
n+1
=(1+h )y
n
设在节点值y
n
上有 δ
n
在y
n+1
上有 δ
n+1
且δ
n+1
仅 δ
n
( 程不的¨ )
~~
n+1 11
(1 ) (1 )
nn
yy hy hyδλλ
++
=?=+?+
~
nn+1
(1 ) (1 ))(
nn
hyy hδ λλ δ∴=+=+?
Euler法?…
+1 nn
||||δ δ?≤ |1 | 1hλ?+ ≤
Euler法的?…的条件 0<h -2/
式Euler?…‰
!式Euler法 y
n+1
=y
n
+hf(x
n+1
,y
n+1
)
于 方程 y’=λy(λ<0)
则y
n+1
=y
n
+hλy
n+1
=> y
n+1
=y
n
/(1-hλ)
δ
δ
λ
n
n+1 n+1 n+1
= y - y =
1-h
~
δ
nnn
=y -y
~
!设y
n
的 则y
n+1
的
! 使?式Euler法?…%|1/(1-hλ)| 1
! λ<0 h>0 上式 成?式Euler法无条件?…
式?…‰
式 y
n+1
=y
n
+h[f(x
n
,y
n
)+f(x
n+1
,y
n+1
)]/2
设 方程 y’=λy(λ<0)
则 y
n+1
=y
n
+h[ y
n
+ y
n+1
]/2
λ
λ
∴
n+1 n
h
1+
2
y= y
h
1-
2
当y
n
有 ′,y
n+1
的
δ
nnn
~
=y -y
n+1 n+1 n+1
y-y =δ
~
λ
δδ
λ
n+1 n
h
1+
2
=
h
1-
2
hλ
1+
2
<1
hλ
1-
2
式?…<=>
λ<0
上式h>0′￠成
式￠?…
HW,
p.118 #16
小结
(())
(())
'
00
y=f x,y
yx =y
Taylor 式
式? ˉ?式ˉ 精ˇ
n+1 n n n
y=y+hf(x,y)
式 一阶
n+1 n n+1 n+1
y=y+hf(x,y)
式 一阶
n+1 n-1 n n
y=y+2hf(x,y) `?式 `阶数值 分法
h
n+1 n n n n+1 n+1
2
y =y + [f(x,y )+f(x,y )]
式 `阶
11
n+1 n 1 222
1nn
2nn1
y=y+k+k
k=hf(x,y)
k=hf(x+h,y+k)
式 `阶
Runge Kutta方法
˙¨?˙¨
收敛与?…‰
Numerical Solutions to Ordinary
Differential Equations
概述
(())
(())
Ι,Problem
00
dy
=f x,y
dx
yx =y
12 1 2 1 2
L,y,y,|f(x,s,y)-f(x,y)| L|y |t,-y ≤
!一阶常微分方程初值问题:
! f(x,y)在[a,b]上连续,且满足Lipschitz条件:
!实际工程技术 生产 科研上会出现大量的微分方程问题很难得到其解析解 有的甚至无法用解析表达式来表示因此只能依赖于数值方法去获得微分方程的数值解求函数y=y(x)满足
y?(x) =f(x,y(x))
则初值问题Problem I有唯一解
!微分方程的数值解法
"不求y=y(x)的精确表达式,而求离散点x
0
,x
1
,…x
n
处的函数值
"设Problem I的解y(x)的存在区间是[a,b] 初始点x
0
=a 取
[a,b]内的一系列节点x
0
,x
1
,...,x
n
a= x
0
< x
1
<…< x
n
=b
一般采用等距
a b
x
0
x
1
x
2
..,x
n-1
x
n
"用数值方法 求得y(x)在 节点x
k
的值y(x
k
) 的
值 用y
k
表示 y
k
y(x
k
) y
0
,y
1
,...,y
n
微分方程的数值解
"求y(x)————>求y
0
,y
1
,...,y
n
! 程
00
n+1 n,n n-1 n-2 n-m
y=y(x)
y=g(h,x y,y,y,...y )
!!方法 采用 式 法
[a,b]n等分 a= x
0
< x
1
<…< x
n
=b h= x
k
=a+kh
b-a
n
10 112 nm nm n n
yy yy yyy
++
→ → →→ → →→ → →
! 式
!Taylor 式
!数值 分法
!,用 微,
≈
n+1 n
nn+1n
y(x )-y(x )
y'(x ) h=x -x
h
00
n+1 n n n
y(x )=y
n=0,1,2,
y=y+hf(x,y)
——? 式 Euler Scheme
4.1? 式
∴≈
n+1 n
n n
y(x )-
y(x
y(x )
f(,)x),
h
(())
(())
'
00
y=f x,y
yx =y
∴≈
n+ n n1 n
y(x )y( +hf(x,) y(xx ))
∴
nnn+1 n
=y +hf(x,yy )
￠￡?￥
1,y(x) 点P
0
(x
0
,y
0
)且在点(x,y)的§currency1
',f(x,y),
2,y(x)在点P
0
(x
0
,y
0
)的§
currency1方程
3,? 2 P
1
f (x
1
,y
1
) '“?'“?y(x)的§currency1 在其上取点
P
2
(x
2
,y
2
) 依此? …
4.ficurrency1P
0
P
1
P
2
P
n
flcurrency1y(x)的
——?ficurrency1法
x
p
0
p
1
p
2
p
3
p
4
x
0
x
1
x
2
x
3
x
4
y
y(x)
(())
(())
'
00
y=f x,y
yx =y
000 0
y=y +f(x,y )(x-x )
在§currency1上取点P
1
(x
1
,y
1
)
01 00
=y + f(xy,yh )
y
1
是Euler 式–求
!用? 微
≈
n+1 n
n+1 n n+1 n+1 n
y(x )-y(x )
y'(x ) y[x,x ]= h=x -x
h
00
n+1 n n+1 n+1
y(x )=y
n=0,1,2,
y=y+hf(x,y)
——?式? 式
法 续
∴≈
n+
n+
1
1n
n+1
y(x )-y(x )
f(x,)y(x )
h
(())
(())
'
00
y=f x,y
yx =y
≈?
n+1 nn1 +n+ 1
+hf(x,y(x ) y() xy(x ))
∴
nnn1 +1+ n+1
=y +hf( yx,y )
here
· 用?式?法 一?
解方程解出y
n+1
的?式表达式,其?…‰
!用? 微
≈
n-1 n+
n+1 n-1
1n
y(x )-y(x )
y'(x ) y[,]=xx
2h
≈
n+1 n-1
n n
y(x )-y(x )
f(x y(,)x)
2h
n-1
00
n+1 nn
y(x )=y
n=0,1,2,
=+2hf(x )y,yy
—— `?法
· ′?用
法求出y
1
用
`?法
法 续
n+1 n n n
y=y+hf(x,y)
式? ˉ?式ˉ
式
n+1 n n+1 n+1
y=y+hf(x,y)
式
n+1 n-1 n n
y=y+2hf(x,y) `?式
˙¨ y(x
n+1
)-y
n+1
˙¨
…￥1?设yn=y(xn),第n 是精确的?,
Rn+1=y(xn+1)-yn+1?法的˙¨,
…法的精ˇ是一阶
…￥2 — 法的˙¨ O(h
p+1
),法有p阶精ˇ.
·
无yn=y(xn)? Rn+1?˙¨
分析 其˙¨ O(h
2
),Taylor 式式分析分析
Euler 式
nn+1 nn
=y +hf(x,yy )
yn=y(xn) y(xn+1)-yn+1 =O(h
2
)
∵
n+1 nnn
=+hf(xy,yy )
nnn
=+hf(y(x ) y,(xx ))
y?(x) =f(x,y(x))
'
=+hy(x) y(x)
T
'
aylor
'
n+1 n n n n+1
2
n
'
y(x )=y(x +h) === y(x )+hy (x )+,(x,x
y)
2
)
(
h ξ
ξ
∈
!
∵
''
2
n+1
() 2
2
n+1
=y( hx)-y h)=O(
y ξ
∴
`?法的˙¨
…￥3?设yn=y(xn),y
n 1
=y(x
n 1
),Rn+1=y(xn+1)-yn+1
`?法的˙¨,
n-1 nn n+1
=+2y hf(x,yy)
…式?法的精ˇ是一阶 `?法的精ˇ是`阶
`?法 其˙¨
yn=y(xn),y
n 1
=y(x
n 1
),
nnn-1
=+2hfy(x ) yx,(x( ))
'
n-1 n
=+y( 2hx) y(x)
Taylor
'
2
''
n+1 n n n n n+
'
1
''
3
n
y(x )=y(x +h)=== y(x )+hy (x )
h
y(x )
2!
++,(
y(
)
)
x,xh
ξ
ξ ∈
3!
∵
n+1 n+
3
1
y(x )-y =O(h )∴
n-1 nn n+1
=+2y hf(x,yy)
'
n-1 n n n
2
n-1
''
3
n n
'
''
(-h)
y(x )
2
y(
y(x )=y(
)
x -h)=y(x )-hy (x )+ +,(x,x )(-h)
!
η
η ∈
3!
上?式 a? ′
''' '''
'
n-1
() ()
3
3
+ nn1
-y(x ) 2hy (x=y(x ) )+ h
yyξη+
!
∴
`?法的˙¨ O(h
3
) 其精ˇ是`阶
,求 的 值
2
(0) 1
dy x
y
dx y
y
=?
=
0.1,0.2,1.0,x = "
解,?
2
(,)
x
fxy y
y
=?
00
0,1,0.1xyh===
式 得:
1
(,)
nn nn
yyhf x y
+
=+
0
1y =
10 00
0
(,)10.1(1 )1.1
1
yyhfxy=+ =+×?=
21 11
20.1
(,) 1.1 0.1 (1.1 ) 1.191818
1.1
yyhfxy
×
=+ =+×? =
32 22
(,) 1.277438yyhfxy=+ =
……
其精确解 21yx=+
¨?o 其?分析
111
()
kkk
e y x y
+++
=?
表? 出¨ 在?
10 9 9 9
( ) (,( )) 1.7330815yyxfxyx=+ =
#
˙¨
10 10
( ) 0.00103yx y?=
#
而¨ 是 10 10
( ) 0.05272yx y?=
数值 分法
(())
(())
'
00
yx
yy
=y
=f x,
n+1 n+1
nn
x
'
x
xx
y( ) f(x,y)xdx dx∴=
∫∫
n+1
n
x
x
1
y(f(x,)(( x)))
nn
yx yx dx
+
=?
∫
a 的… 分用数值 分 式求 值
1 用 数值 分 式
n+1
n
1
x
x
(f(x,) (,)y(x) y() )
nnn n
dfxx xxx
+
≈
∫
1
y( )-y( ) ( y(,))
nn nn
hf xxx x
+
≈∴
1
y( )y( ) + (,)y( )
nnnn
xhfxx
+
≈
+(,)
1
yyy
n
hf x
n nn
↓↓ ↓
=
+
2 用 式
1
1
n+
n
nn
x
nn+11
(
+
)
2
x
[y(xf(x,) f(x,))yy(x)+f(x,)]
nn
xx
dx
+
≈
∫
n+1 nn+
2
1nn+1n
yf(x,) fy(x )-y (x,(x ) y(x(x ) [ )])+
h
∴≈
nnn +1n+ +
2
1 1
y=y [f(x +],) f(x,)yy
h
∴+
—— 式
n+1 n n+1 n+1
y=y+hf(x,y)
n+1 n n n
y=y+hf(x,y)
# 式,?示? 式,?式? 式? 得
# 式是?式? 式,其精ˇ `阶
式的精ˇ
y
n
=y(x
n
),Talor 式有分析 式是?式 式 其˙¨ O(h
3
)
用 用到到`?函数的Taylor 式
…? 式 的精ˇ是2阶的.
n nn+1n1
2
+1+
yf=y [ (x,) f(x,)y+]y
h
+
f(x
n+1
,y
n+1
)=f(x
n+1
,y(x
n+1
)+(y
n+1
-y(x
n+1
))
=f(x
n+1
,y(x
n+1
))+f
y
(x
n+1
,)(y
n+1
-y(x
n+1
)) (x
n
x
n+1
)
=y’(x
n+1
)+f
y
(x
n+1
,)(y
n+1
-y(x
n+1
))
=y’(x
n
)+hy”(x
n
)+O(h
2
) +f
y
(x
n+1
,)(y
n+1
-y(x
n+1
))
=f(x
n
,y
n
)+hy”(x
n
) +f
y
(x
n+1
,)(y
n+1
-y(x
n+1
)) +O(h
2
)
y(x
n+1
)= y(x
n
+h) = y(x
n
)+hy’(x
n
)+h
2
y”(x
n
) /2 +O(h
3
)
=y
n
+hf(x
n
,y
n
)+h
2
y”(x
n
)/2+O(h
3
)
=y
n
+hf(x
n
,y
n
)/2+h[f(x
n
,y
n
)+hy”(x
n
)]/2+O(h
3
)
而y(x
n+1
)=y
n+1
+h f
y
(x
n+1
,)[y(x
n+1
)-y
n+1
] /2 +O(h
3
)
y(x
n+1
)-y
n+1
=h f
y
(x
n+1
,)[y(x
n+1
)-y
n+1
] /2 +O(h
3
)
y(x
n+1
)-y
n+1
=O(h
3
)/[1-hf
y
(x
n+1
,)/2]=O(h
3
)
式的?˙¨ O(h
3
) 其精ˇ是精ˇ是2阶
f(x
n+1
,y
n+1
)=f(x
n
,y
n
)+hy”(x
n
) +f
y
(x
n+1
,)(y
n+1
-y(x
n+1
)) +O(h
2
)
y(x
n+1
)= y
n
+hf(x
n
,y
n
)/2+h[f(x
n
,y
n
)+hy”(x
n
)]/2+O(h
3
)
= y
n
+hf(x
n
,y
n
)/2+h [f(x
n+1
,y
n+1
)
-f
y
(x
n+1
,)(y
n+1
-y(x
n+1
)) +O(h
2
)]/2+O(h
3
)
= y
n
+h[f(x
n
,y
n
) + f(x
n+1
,y
n+1
) ]/2
+ h f
y
(x
n+1
,)(y(x
n+1
) - y
n+1
) /2 +O(h
3
)
式的?用
dy
=y,y(0)=1.
dx
4.1 用 式求初值问题的解在x=0.01上的值y(0.01).
10 01
h
y=y+ [y+y]
2
解,取h=0.01 x
0
=0 y
0
=y(0)=1 则 y(0.01) y
1
f(x,y)=y 式,
≈+ + =
23
0.01
0.01 0.01
y(0.01)=e =1+0.01+ + +...
2! 3!
.
..
!
2
001
1 0 01 1 01005
2
n+1 n n n n+1 n+1 n n n+1
hh
y =y + [f(x,y )+f(x,y )]=y + [y +y ]
22
h
2
h
0
2
1
1
=
-
1+
yy
≈
10
2
2
2
h
y=(1+ )( )y
2
hh
(1+ ) + =1.01005
2
hh
1+ + +...
24
4
于?…‰
解析解
x
y=e
式的?
精ˇ有–,”?式,? 求解,
量大
式法
式,” 初值,且第2 初值只能 其
方法给出,能?o的 精ˇ有影响
`?法
…‰?式?法最简?,精ˇ低
法
4.2 改 的的Euler法
!Euler 式 y
n+1
=y
n
+hf(x
n,
y
n
)?式 一阶
! 式 y
n+1
=y
n
+h[f(x
n
,y
n
)+f(x
n+1
,y
n+1
)]/2?式 `阶
!Euler 式 量小 精ˇ低
! 式 量大 ˇ
综合? 式 出预报—校 式
n
_
n+1 nn
=y +hf(x,yy )
预报
h
nnnn
_
n
n+12
+1
+1
=y + [f(x,y )+f(x,]y )
y
校
——改 的Euler法写成嵌套 式 n+1
y
h
nnnn+12 nnn
=y + [f(x,y )+f(x y+hf(x,y,))]
式? 法
化 式
n+1
y
nnn nnn+1 n
1
2 n
y+y+hf(x,y) y+hf(=[ x,yhf(x,+ ))]
pn nn
cn n+1p
1
n+1 p c
2
y=y+hf(x+y)
y=y+hf(x +y)
y=(y+y)
4.4 用改 的Euler法解初值问题在区间[0,0.4]上
h=0.1的解 并? 与精确解的 异
2
dx
=y
dy
y(0) = 1
说 精确解y=1/(1-x)
解 Euler法的具 式 y
n+1
=y
n
+hy
n
2
改 的Euler法的具 式
2
pn n
2
cn p
1
n+1 p c
2
y=y+hy
y=y+hy
y=(y+y)
x
0
=0,h=0.1,则
x
1
=0.1,x
2
=0.2,x
3
=0.3,x
4
=0.4
y
1:
y
p
=y
0
+0.1y
0
2
=1+0.1·1
2
=1.1
y
c
=1+0.1X1.1
2
=1.121
y
1
=(1.1+1.121)/2 1.1118
求y
2,
y
3,
y
4
见P95表:
y(x
n
)-y
n
随n 大而 大 表 ¨?
·
1 y(x
n
)=y
n
导改 的Euler法的˙¨
O(h
3
),具有`阶精ˇ
2 改 的Euler法也 写成如 化 式
n+1 n 1 2
1nn
2n+1,n1
1
y=y+(k+k)
2
k=hf(x,y)
k=hf(x y+k)
n+1
y
h
nnnn+12 nnn
=y + [f(x,y )+f(x y+hf(x,y,))]
HW,
p.117 #3 4
龙格—库塔方法
1n+1 n
1nn
y=y+
k=hf(x,y)
k
改 的Euler 式,
1
12n+1 n
1nn
2n+1
2
,n 1
y=y+
k=hf(x,y)
k=hf(x y+k)
(k +k )
'(,)
()
00
yfxy
yx y
=
=
Euler 式 y
n+1
=y
n
+hf(x
n,
y
n
)
写成精ˇ 一阶精ˇ ``阶阶
n+1 n
(x ) (x )
y-y
h
Lagrange?值…?
ξ?∈
nn+1
(x,x )
ξ
'
=y( )
ξ∴
n+1 n
'
y(x ) = y( + hy) (x )
ξξξ而
'
=f((),y(y ))
n+1
*
n
y(x ) = (x ) +kyξξ∴
n+1 n
hf(y(x ) = y,y()+ (x ))
y(x)在[x
n
,x
n+1
]上的? '“
*
nn 1
k=hf(x,y)=k
!取
—Euler 式 !取
12
kk
*
22
k= +
—改 Euler 式
Euler 式用一 点的值k
1
k
*
的 值 而改 的
Euler 式用` 点的值k
1
k
2
的? 值? k
*
值 改
的Euler法?Euler法精ˇ
龙格 库塔法
Runge-Kutta法的 在 x
n,
x
n+1
内多预报￠ 点的k
i
值并用其 权? 值? k
*
值 而构造出具有更 精ˇ的
式
4.3.2 `阶龙格—库塔方法
k
1
与k
2
的 权? 来 k
*
设
αβ
n+1 n 1 2
1n
12
n
2n n1
y=y+k+k
k = hf(x,y ) (*)
k=hf(x+h,y+ k)
ww
其?w
1,
w
2,
α,β 待…参数 适当选取参数 使(*)式的精ˇ `阶 使其˙¨ O(h
3
)
y(x
n
)=y
n
泰勒 式
2
''' 3
h
n+1n n n2
y(x )=y(x )+hy(x )+ y (x )+O(h )
∵
'
nnn nn
y (x ) = f(x,y(x )) = f(x,y )
'' ' ' ' '
xxy xy
( y (x)=(y (x)) =[f(x,y)] =f +f y (x)=f +f f(x,y) )∵
''
nxnnynn nn
y (x )=f (x,y )+f (x,y )f(x,y )
∴
2
3
h
nnn xnn+1 n ynn nn2
=y+hf(x,y)+ [f(x,y)y(x +f (x,y )f(x,y )]+O(h) )
`阶龙格 库塔法
αβα β
2
2 n 1 nn xnn 1ynn
k =hf(x + h,yn+ k )=h{f(x,y )+ hf (x,y )+ k f (x,y )+O(h )}
多?函数的泰勒 式
αβ
n+1 n 1 2
1n
12
n
2n n1
y=y+k+k
k = hf(x,y ) (*)
k=hf(x+h,y+ k)
ww
α
β
∴
2
n1nn 2nn 2xnn
23
2ynnn
+1
n
n
=y +w hf(x,y )+w hf(x,y )+w h f (x,y )
+w h f (x,y )f(x,y )+O(h
y
)
α
β
2
h
nnxn2n+ 2
3
ynn nn2
1 12
=y+ hf(x,y)+ [ f(x,y)
+ f (x,y )f(x,y )]+O(h )
2
(w +w )
w
2w
y
(2)
2
3
h
nnn xnnynnnn2n+1
=y +hf(x,y )+ [f (x,y )+f (x,y )f(x,y )]+O(h ) y(x ) (1)
则有
w
1
+w
2
=1
w
2
α=1/2
w
2
=1/2
· 上述方程组有四 未知量 只有三 方程 有无穷多组解
<1>与<2> 使 y(x
n+1
)-y
n+1
=O(h
3
)
`阶龙格 库塔法
!取一组解便得一种`阶龙 库 式
αβ
n+1 n 1 2
1nn
2
12
nn1
y=y+k+k
k=hf(x,y)
k=hf(x
w
+h,y+k)
w
w
1
+w
2
=1
w
2
α=1/2
w
2
=1/2
!当w
1
=w
2
=1/2,= =1′`阶Runge-Kutta 式
y
n+1
=y
n
+k
1
/2+k
2
/2
k
1
=hf(x
n
,y
n
)
k
2
=hf(x
n
+h,y
n
+k
1
)
此 改 得Euler法
!取w
2
=0 w
2
=1 =1/2 =1/2
y
n+1
=y
n
+k
2
k
1
=hf(x
n
,y
n
)
k
2
=hf(x
n
+h/2,y
n
+k
1
/2)
此?点法?变 的 Euler 式三阶龙格 库塔法
!三阶龙格 库塔法是用三 值k
1
,k
2
,k
3
的 权? 来 k*
有
y
n+1
=y
n
+c
1
k
1
+c
2
k
2
+c
3
k
3
k
1
=hf(x
n
,y
n
)
k
2
=hf(x
n
+a
2
h,y
n
+b
21
k
1
)
k
3
=hf(x
n
+a
3
h,y
n
+b
31
k
1
+b
32
k
2
)
! 使其具有三阶精ˇ 必须使˙¨ O(h
4
)
!? `阶龙格 库塔法的 导 c
1
,c
2
,c
3
,a
2
,a
3
,b
21
,b
31
,b
32
满足
c
1
+c
2
+c
3
=1
a
2
=b
21
a
3
=b
31
+b
32
c
2
a
2
+c
3
a
3
=1/2
c
2
a
2
2
+c
3
a
3
2
=1/3
c
3
a
3
2
=1/6
方程组解
得三阶龙格-库塔 式
Kutta 式
k
n+1
=y
n
+(k
1
+4k
2
+k
3
)/6
k
1
=hf(x
n
,y
n
)
k
2
=hf(x
n
+h/2,y
n
+k
1
/2)
k
3
=hf(x
n
+h,y
n
-k
1
+2k
2
)
四阶龙格 库塔法
!? 出四阶龙格-库塔 式 常用的有
经典Runge-Kutta法
y
n+1
=y
n
+(k
1
+2k
2
+2k
3
+k
4
)/6
k
1
=hf(x
n
,y
n
)
k
2
=hf(x
n
+h/2,y
n
+k
1
/2)
k
3
=hf(x
n
+h/2,y
n
+k
2
/2)
k
4
=hf(x
n
+h,y
n
+k
3
)
˙¨ O(h
5
)
!还有,Gill 式及m (m>4)阶龙格 库塔法
m>4′ 量太大 精确ˇ不一… 有′会降低龙格 库塔法
于经典的四阶Runge-Kutta法给出如? 法求解 dy/dx=f(x,y) a x b
y (a)=y
0
!Step 1,输入a,b,y
0
及N
!Step 2,(b-a)/N=>h,a=>x,y
0
=>y
!Step 3,输出 (x,y)
!Step 4,For I=1 T0 N
" hf(x,y)=>k
1
" hf(x+h/2,y+ k
1
/2)=> k
2
" hf(x+h/2,y+k
2
/2)=>k
3
" hf(x+h,y k
3
)=>k
4
" y+(k
1
+2k
2
+2k
3
+k
4
)/6=>y
" x+h=>x
" 输出(x,y)
!END
龙格 库塔法
用四阶经典Runge Kutta方法解初值问题
()
=
=
10
2
y
y
x
y
dx
dy
2.0=h
(1)求
1
y
2.0,1,0
00
=== hyx
() 1
2
,
0
0
0001
=?==
y
x
yyxfK
9181818.0
2
2
2
2
2
,
2
10
0
10
1002
=
+
+
+=
++=
K
h
y
h
x
K
h
y
K
h
y
h
xfK
9086375.0
2
2
2
2
2
,
2
20
0
20
2003
=
+
+
+=
++=
K
h
y
h
x
K
h
y
K
h
y
h
xfK
()
()
8432399.0
2
,
30
0
30
3004
=
+
+
+=
++=
hKy
hx
hKy
hKyhxfK
()
1832293.1
2
6
432101
=
++++= KKKK
h
yy
() 1832160.14.112
11
==+= xxy
5
1
103.1
×=e
龙格 库塔法
(2)求
2
y
1
0.2,0.2xh==
1832293.1
1
=y
() 8451714.0
2
,
1
1
1111
=?==
y
x
yyxfK
()
2111
22
,0.7946656
hh
Kfx yK=++=
()
3112
22
,0.7874945
hh
Kfx yK=++=
()
4113
,0.7440375KfxhyhK=++=
()
21 1 2 3 4
2 1.3416803
6
h
yy KK KK=+ ++ + =
() 3416408.112
22
=+= xxy
5
2
100.4
×=e
变 龙格 库塔法
!用户 出问题I
问题 如￡判˙|y(x
n
)-y
n
|<ε
精确值y(x
n
)未知如￡取h=?
解 如用p阶龙格 库塔法˙¨ O(h
p+1
)
y’=f(x,y)
y(x
0
)=y
0
求¨ < 10
-8
求数值解
x
n
h/2 x
n+1
h
!如 x
n
-----------------$x
n+1
y
n
=y(x
n
)y
n+1
(h)
则 y(x
n+1
)-y
n+1
(h)
ch
p+1
fi半x
n
$x
n+h/2
$x
n+1
分? y(x
n+1
)的 值y
n+1
(h/2)
则 y(x
n+1
)-y
n+1
(h/2)
2c(h/2)
p+1
h
2
n+1
()
n+1
(h) p
n+1 n+1
y(x )-y
1
y(x )-y 2
∴≈
h
2
h
2
nn +1+1
( ()
(h)
)
n+1
p
n +1
1
[]
2
y-yyy
1
(x )- ≈
$
变 龙格 库塔法
…? 于问题I—用P阶龙格 库塔法 y(x
n+1
)在 fi半
的 值分别 y
n+1
(h)
,y
n+1
(h/2)
则有¨ 式
p
h
2
n+1
h
2
n+1
1
2
()
n
()
(h)
n+
1
1 1+
y(x || )-y yy|- =|
≈?
· 1
0
¨ 的事?估 法
2
0
停机准则 <ε ( 保 |y(x
n+1
)-y
n+1
(h/2)
|< )
解
h取大˙¨ ch
p+1
大 不精确
h取小 运 量大 多 舍入¨?大解决策略 变 龙格 库塔法
If > fi半反复 直至 < 止取h 最?一次的 y
n+1
最?一次 的结果
Else if ( < ) 倍反复 直至 > 止最?一次运 的 一次 结果 –?
HW,p.117 #4
用经典四阶Runge-Kuttta法
4.5 收敛与?…‰
于常微分方程初值问题
!求y=y(x)--------$求y(x
n
)------------$求y
n
x
n
=x
0
+nh
离散化 种 式
Euler法 改 的Euler法 龙格 库塔法 是? 法数值解法
(())
(())
'
00
y=f x,y
yx =y
法 y
n+1
′只用到 一 的结果y
n
!?式? 法
y
n+1
=y
n
+hφ(x
n
,y
n
,h)
(x,y,h) 量函数 它依赖于f,仅是x
n
,y
n
,h的函数
!Def:— 数值方法 于固…的x
n
=x
0
+nh,当
h->0(n->∞)′ y
n
->y(x
n
) 则 法是收敛的
0
(0) )(
lim
n
h
n
n
yx y
→
→∞
=
(x
n
=x
0
+nh 固…值)
收敛与?…‰
· 数值方法是ˉ收敛取决于¨ y(x
n
)-y
n
的变化情况
于p阶 式,其˙¨ O(h
p+1
)
y(x
n+1
)-y
n+1
=O(h
p+1
) 其?…y
n
=y(x
n
).
虽h->0′˙¨ y(x
n+1
)-y
n+1
->0 并不能说
法收敛 因其 y
n
=y(x
n
)不成 此 入——
Def,y(x
n
)-y
n
法的解y
n
的?˙¨
˙¨?˙¨ 的区别是
法收敛%
0
(0) )(
lim
n
h
n
n
yx y
→
→∞
=
h->0′?˙¨ y(x
n+1
)-y
n+1
->0
Φ
ΦΦ≤
__ __
|(x,y,h)- (x,y,h)|L|y-y|
Th —? 法y
n+1
y
n
+h φ(x
n
,y
n
,h)具有p阶精ˇ 且 量函数
φ(x,y,h) 于y满足Lipschitz条件
—初值y
0
是准确的(y
0
=y(x
0
)) 则其?˙¨
y(x
n
)-y
n
=O(h
p
).
·
—? 法满足?上条件,? 其收敛改 Euler法的收敛‰
!判…改 Euler法的收敛‰
y
n+1
=y
n
+h[f(x
n
,y
n
)+f(x
n
+h,y
n
+hf(x
n
,y
n
))]/2
则 (x,y,h)=[f(x,y)+f(x+h,y+hf(x,y))]/2
— y
o
=y(x
0
)
≤
__ __
|f(x,y)-f(x,y)| L|y-y|
f(x,y) 于y满足 --条件:
1
2
ΦΦ
≤
__ __ _
__
_
|(x,y,h)- (x,y,h)|
[|f(x,y)- f(x,y)|+|f(x+h,)- f(x+y+hf(x,y) y+hh,f(x,y))|]
则
1
2
≤
__ __ __
[L|y- y |+L|y+hf(x,y)- y-hf(x,y)|]
≤
__ __ __ __
2
1
[L|y- y |+L|y- y |+L h
h
L(1+ L)
2
|y-y|]= |y-y|
2
…h h
0
则
ΦΦ≤
__
0
__
|(x,y,h)- (x,y,h
h
L(1+ L)
2
)| |y-y|
(x,y,h) 满足 条件,Th改 的Euler法收敛
! —f(x,y) 于y满足 条件,且 y
0
=y(x
0
)′
Euler法 准四阶龙格 库法也收敛
…‰
! 在 收敛‰′ 一般 数值方法 程是准确的 实际上 并 如此初始值y
0
有¨ δ
0
y
0
-y(x
0
).
的 一 有舍入¨
! 初始¨ 在 程的? 是 的 还是
‰ 是?…‰问题
Def 设在节点x
n
处用数值法得到的? 数值解 y
n
而实际
得到的 值 第 第n 数值解的
n
~
y
δ
nnn
~
=y -y
Def —一种数值方法在节点x
n
处的数值解y
n
的 δ
n
0
而在的 节点值y
m
(m>n)上有 δ
m
.
当|δ
m
| |δ
n
| (m=n+1,n+2,……) 则 数值 法是?…的
! 分析 法的数值?…‰ 常 方程
y’=λy (λ<0)
Euler法的?…‰
Euler法 y
n+1
=y
n
+hf(x
n
,y
n
)
方程 y’=λy (λ<0)
y
n+1
=(1+h )y
n
设在节点值y
n
上有 δ
n
在y
n+1
上有 δ
n+1
且δ
n+1
仅 δ
n
( 程不的¨ )
~~
n+1 11
(1 ) (1 )
nn
yy hy hyδλλ
++
=?=+?+
~
nn+1
(1 ) (1 ))(
nn
hyy hδ λλ δ∴=+=+?
Euler法?…
+1 nn
||||δ δ?≤ |1 | 1hλ?+ ≤
Euler法的?…的条件 0<h -2/
式Euler?…‰
!式Euler法 y
n+1
=y
n
+hf(x
n+1
,y
n+1
)
于 方程 y’=λy(λ<0)
则y
n+1
=y
n
+hλy
n+1
=> y
n+1
=y
n
/(1-hλ)
δ
δ
λ
n
n+1 n+1 n+1
= y - y =
1-h
~
δ
nnn
=y -y
~
!设y
n
的 则y
n+1
的
! 使?式Euler法?…%|1/(1-hλ)| 1
! λ<0 h>0 上式 成?式Euler法无条件?…
式?…‰
式 y
n+1
=y
n
+h[f(x
n
,y
n
)+f(x
n+1
,y
n+1
)]/2
设 方程 y’=λy(λ<0)
则 y
n+1
=y
n
+h[ y
n
+ y
n+1
]/2
λ
λ
∴
n+1 n
h
1+
2
y= y
h
1-
2
当y
n
有 ′,y
n+1
的
δ
nnn
~
=y -y
n+1 n+1 n+1
y-y =δ
~
λ
δδ
λ
n+1 n
h
1+
2
=
h
1-
2
hλ
1+
2
<1
hλ
1-
2
式?…<=>
λ<0
上式h>0′￠成
式￠?…
HW,
p.118 #16
小结
(())
(())
'
00
y=f x,y
yx =y
Taylor 式
式? ˉ?式ˉ 精ˇ
n+1 n n n
y=y+hf(x,y)
式 一阶
n+1 n n+1 n+1
y=y+hf(x,y)
式 一阶
n+1 n-1 n n
y=y+2hf(x,y) `?式 `阶数值 分法
h
n+1 n n n n+1 n+1
2
y =y + [f(x,y )+f(x,y )]
式 `阶
11
n+1 n 1 222
1nn
2nn1
y=y+k+k
k=hf(x,y)
k=hf(x+h,y+k)
式 `阶
Runge Kutta方法
˙¨?˙¨
收敛与?…‰
