Download Analysis of Numerical Methods for Solving Ordinary Differential Equations and more Study notes Mathematics in PDF only on Docsity!
- Thurs 1/6 1
- IVP for ODEs y ′^ = f (y) , y(0) = y 0 : find y(t) ex y ′^ = y , y(0) = y 0 ⇒ y(t) = y 0 et y ′^ = sin y , y(0) = y 0 ⇒ y(t) =? def An IVP is well-posed if the following conditions are satisfied.
- A solution exists.
- The solution is unique.
- The solution depends continuously on the data (i.e. y 0 , f ). def f satisfies a Lipschitz condition on a domain D if there exists a constant L st |f (y) − f (u)| ≤ L |y − u| for all y , u ∈ D. note : we typically assume that f (y) is smooth, so L = max |fy | thm If f satisfies a Lipschitz condition, then the IVP is well-posed. ex
y ′^ = ay + b , y(0) = y 0 ⇒ L = |a| , y(t) = y 0 eat^ + ba (eat^ − 1)
y ′^ = y 1 /^2 , y(0) = 0 ⇒ L = ∞ for y ≥ 0 , non-unique solution , y(t) = 0 , t
2 4
y ′^ = y 2 , y(0) = 1 ⇒ L = ∞ for y ≥ 1 , y(t) = (^1) −^1 t , lim t→ 1 y(t) = ∞
Euler’s method h = ∆t = step size , tn = nh y (^) n = y(tn ) : exact solution , u (^) n : numerical approximation y ′^ = f (y) : differential equation u (^) n+1 − u (^) n h =^ f^ (u^ n^ )^ :^ finite-difference scheme u (^) n+1 = u (^) n + hf (u (^) n ) , u 0 = y (^0) ex y ′^ = y , y 0 = 1 u 0 = 1 u 1 = u 0 + hf (u 0 ) = 1 + h u 2 = u 1 + hf (u 1 ) = (1 + h) + h(1 + h) = (1 + h)^2 · · · u (^) n = (1 + h)n tn = 1 ⇒ nh = 1 , y (^) n = y(1) = e = 2. 7182818...
n h u (^) n y (^) n − u (^) n (y (^) n − u (^) n )/h 10 0.1 2.5937425 0.1245 1. 20 0.05 2.6532977 0.0650 1. 40 0.025 2.6850638 0.0332 1. 80 0.0125 2.7014849 0.0168 1.
convergence proof (special case) y ′^ = y , y 0 = 1 ⇒ y(t) = et u (^) n+1 = u (^) n + hun , u 0 = 1 ⇒ u (^) n = (1 + h)n
consider t = 1 , h = n^1
h^ lim→ 0 u^ n^ =^ nlim→∞
( 1 +^1 n
)n = (^) nlim→∞ exp
( n ln
( 1 + n^1
)) = exp
( n^ lim→∞ n^ ln
( 1 + n^1
))
n^ lim→∞ n^ ln
( 1 + n^1
) = ∞ · 0 = (^) nlim→∞^ ln (1 +^ n^1 ) n^1 =^ nlim→∞
1 + (^1) n^ ·^
− n 21 · (^) −^11 n 2
⇒ lim h→ 0 u (^) n = e = y(1) ⇒ Euler’s method converges
hw : y ′^ = ay + b , y(0) = y (^0) convergence proof (general case) y ′^ = f (y) , y(0) = y (^0) u (^) n+1 = u (^) n + hf (u (^) n ) , u 0 : given fix t > 0 , set h = t/n ⇒ t = nh = tn , y (^) n = y(tn ) goal : (^) hlim→ 0 u (^) n = y (^) n
y (^) n+1 = y (^) n + hf (y (^) n ) + τ (^) n , τ (^) n : local truncation error step 1 : estimate τ (^) n y (^) n+1 = y(tn+1 ) = y(tn + h) = y(tn ) + hy ′^ (tn ) + h
2 2 y^
′′ (^) ( ˜t )
y (^) n+1 = y (^) n + hf (y (^) n ) + h
2 2 y^
′′ (^) ( ˜t ) ⇒ τ (^) n = h^2 2 y^
′′ (^) ( ˜t )
step 2 : analyze error , en = y (^) n − u (^) n : error en+1 = y (^) n+1 − u (^) n+1 = y (^) n + hf (y (^) n ) + τ (^) n − (u (^) n + hf (u (^) n )) = y (^) n − u (^) n + h(f (y (^) n ) − f (u (^) n )) + τ (^) n = en + h(f (y (^) n ) − f (u (^) n )) + τ (^) n
|en+1 | ≤ |en | + h|f (y (^) n ) − f (u (^) n )| + |τ (^) n | ≤ |en | + hL|y (^) n − u (^) n | + τ , τ = max |τ (^) n | |en+1 | ≤ (1 + hL)|en | + τ |e 1 | ≤ (1 + hL)|e 0 | + τ |e 2 | ≤ (1 + hL)|e 1 | + τ ≤ (1 + hL)^2 |e 0 | + (1 + (1 + hL))τ |e 3 | ≤ (1 + hL)|e 2 | + τ ≤ (1 + hL)^3 |e 0 | + (1 + (1 + hL) + (1 + hL)^2 )τ · · · |en | ≤ (1 + hL)n^ |e 0 | + n j∑−=0^1 (1 + hL)j^ τ
n∑− 1 j=0^ (1 +^ hL)
j (^) = (1 +^ hL)n^ −^1 (1 + hL) − 1 =^
(1 + hL)n^ − 1 hL note : 1 + x ≤ ex^ ⇒ (1 + x)n^ ≤ enx^ ⇒ (1 + hL)n^ ≤ enhL^ = eLt
|en | ≤ eLt^ |e 0 | + e Lt (^) − 1 hL τ
τ = max |τ (^) n | = h
2 2 M^ ,^ M^ = max^ |y^
′′ (^) (t)|
|en | ≤ eLt^ |e 0 | + e
Lt (^) − 1 L
M h 2 if lim h→ 0 |e 0 | = 0 , then lim h→ 0 |en | = 0 ⇒ (^) hlim→ 0 u (^) n = y (^) n , i.e. Euler’s method converges note There are two types of convergence.
- Fix t, set h = t/n. Does lim h→ 0 u (^) n = y(t)?
- Fix h. Does limn→∞ u (^) n = y(∞)?
ex y ′^ = y , y(0) = 1 E ′^ = E − 12 et^ ⇒ E(t) = − 12 tet^ , E(1) = − 1. 36 ⇒ u (^) n − y (^) n = − 1. 36 h + O(h^2 ) pf define dn = u (^) n − (y (^) n + hEn ) we already know that dn = O(h) , we must show that dn = O(h^2 ) dn+1 = u (^) n+1 − (y (^) n+1 + hEn+1 ) = u (^) n + hf (u (^) n ) − (y (^) n + hy ′ n + 12 h^2 y (^) n′′ + O(h^3 ) + h(En + hE (^) n′ + O(h^2 ))) = u (^) n − (y (^) n + hEn ) + h(f (u (^) n ) − y ′ n ) − h^2 ( 12 y ′′ n + E (^) n′ ) + O(h^3 ) f (u (^) n ) = f (y (^) n + hEn + dn ) = f (y (^) n ) + fy (y (^) n )(hEn + dn ) + O(h^2 ) f (u (^) n ) − y (^) n′ = fy (y (^) n )(hEn + dn ) + O(h^2 ) dn+1 = dn + h(fy (y (^) n )(hEn + dn ) + O(h^2 )) − h^2 ( 12 y (^) n′′ + E (^) n′ ) + O(h^3 ) = dn + hfy (y (^) n )dn + h^2 (fy (y (^) n )En − 12 y (^) n′′ − E (^) n′ ) + O(h^3 ) ︸ ︷︷ ︸ (^) = 0 by definition of E(t)
dn+1 = dn + hfy (y (^) n )dn + O(h^3 ) |dn+1 | ≤ (1 + hL)|dn | + κ , κ = O(h^3 ) , d 0 = 0 · · · |dn | ≤ e
Lt (^) − 1 hL κ^ =^ O(h
(^2) ) ok
note u (^) n − y (^) n = O(h) u (^) n − y (^) n = hEn + O(h^2 ) u (^) n − y (^) n = hEn + h^2 Dn + O(h^2 ) , where Dn = D(tn ) hw : find the equation satisfied by D(t)
application : Richardson extrapolation u (^) n = y (^) n + hEn + h^2 D (^) n + · · · A 00 = u h^ = y + hE + h^2 D + O(h^3 ) A 10 = u h/^2 = y + h 2 E + ( h 2 )^2 D + O(h^3 )
A 20 = u h/^4 = y + h 4 E + ( h 4 )^2 D + O(h^3 ) eliminate O(h) term 2 A 10 − A 00 = A 11 = y − h 2 2 D + O(h^3 )
2 A 20 − A 10 = A 21 = y − h 8 2 D + O(h^3 ) eliminate O(h^2 ) term 4 A 21 − A (^11) 3 =^ A^22 =^ y^ +^ O(h
note A (^) j 0 = y + O(h)
A (^) j 1 = 2A (^) j 0 − A (^) j− 1 , 0 = y + O(h^2 )
A (^) j 2 = 4 A^ j^1 − 3 A^ j−^1 ,^1 = y + O(h^3 )
A (^) j 3 = 8 A^ j^2 − 7 A^ j−^1 ,^2 = y + O(h^4 )
ex : y ′^ = y , y(0) = 1 ⇒ y(1) = e = 2. 7182818 j h A (^) j 0 A (^) j 1 A (^) j 2 A (^) j 3 0 0.1 2. 1 0.05 2.6532977 2. 2 0.025 2.6850638 2.7168299 2. 3 0.0125 2.7014849 2.7179060 2.7182647 2. O(h) O(h^2 ) O(h^3 ) O(h^4 ) down a column : decreasing step size , fixed order of accuracy , h-refinement across a row : fixed step size , increasing order of accuracy , p-refinement
Runge-Kutta methods y ′^ = f (y) k 1 = f (u (^) n ) k 2 = f (u (^) n + ahk 1 ) u (^) n+1 = u (^) n + h(bk 1 + ck 2 ) : 2-stage RK method choose a , b , c to minimize the local truncation error y (^) n+1 = y (^) n + h(bf (y (^) n ) + cf (y (^) n + ahf (y (^) n ))) + τ (^) n y (^) n+1 = y (^) n + hy (^) n′ + h 2 2 y (^) n′′ + O(h^3 ) f (y (^) n ) = y (^) n′ f (y (^) n + ahf (y (^) n )) = f (y (^) n ) + fy (y (^) n ) · ahf (y (^) n ) + O(h^2 ) = y (^) n′ + ahy (^) n′′ + O(h^2 ) h(by (^) n′ + c(y ′ n + ahy ′′ n + O(h^2 )) + τ (^) n = hy ′ n + h 2 2 y ′′ n + O(h^3 ) h(b + c − 1)y ′ n + h^2 (ac − 12 )y (^) n′′ + τ (^) n = O(h^3 ) b + c − 1 = 0 , ac − 12 = 0 ⇒ τ (^) n = O(h^3 )
c = (^21) a , b = 1 − (^21) a : 1-parameter family of methods
midpoint method : a = 12 , b = 0 , c = 1 k 1 = f (u (^) n ) k 2 = f (u (^) n + h 2 k 1 ) u (^) n+1 = u (^) n + hk 2 = u (^) n + hf (u (^) n + h 2 f (u (^) n )) modified Euler method : a = 1 , b = 12 , c = (^12) k 1 = f (u (^) n ) k 2 = f (u (^) n + hk 1 ) u (^) n+1 = u (^) n + h 2 (k 1 + k 2 ) = u (^) n + h 2 (f (u (^) n ) + f (u (^) n + hf (u (^) n )))
4th order Runge-Kutta k 1 = f (u (^) n ) k 2 = f (u (^) n + h 2 k 1 ) k 3 = f (u (^) n + h 2 k 2 ) k 4 = f (u (^) n + hk 3 ) u (^) n+1 = u (^) n + h 6 (k 1 + 2k 2 + 2k 3 + k 4 ) : 4 stages ex y ′^ = y , y(0) = 1 k 1 = u (^) n k 2 = u (^) n + h 2 k 1 = (1 + h 2 )u (^) n k 3 = u (^) n + h 2 k 2 = (1 + h 2 (1 + h 2 ))u (^) n = (1 + h 2 + h 4 2 )u (^) n
k 4 = u (^) n + hk 3 = (1 + h(1 + h 2 + h 4 2 ))u (^) n = (1 + h + h 2 2 + h 4 3 )u (^) n
u (^) n+1 = u (^) n + h 6 (1 + 2(1 + h 2 ) + 2(1 + h 2 + h 4 2 ) + (1 + h + h 2 2 + h 4 3 ))u (^) n u (^) n+1 = (1 + h + h 2 2 + h 6 3 + 24 h^4 )u (^) n t = 1 , h = 1/n h u (^) n y (^) n − u (^) n (y (^) n − u (^) n )/h^4 0.2 2.7182511 0.00003069185 0. 0.1 2.7182797 0.00000208432 0. 0.05 2.7182817 0.00000013580 0.
- Tues 1/18 13
note
- (a) ⇒ the difference scheme approximates the correct differential equation u (^) n+1 = u (^) n + hF (u (^) n , h) y (^) n+1 = y (^) n + hF (y (^) n , h) + τ (^) n y (^) n+1 − y (^) n h =^ F^ (y^ n^ , h) +^
τ (^) n h ⇒^ y^ n′ =^ f^ (y^ n^ ) In this case we say that the difference scheme is consistent with the differential equation.
- The scheme is stable wrt initial data if |u (^) n − v (^) n | ≤ C|u 0 − v 0 |, where u (^) n , v (^) n are the numerical solutions starting from u 0 , v 0 and C is independent of h, i.e. a small change in the initial data of the difference scheme leads to a small change in the solution of the difference scheme.
- The theorem says that consistency + stability ⇒ convergence. pf
- u (^) n+1 = u (^) n + hF (u (^) n , h) y (^) n+1 = y (^) n + hF (y (^) n , h) + τ (^) n y (^) n + hy ′ n + O(h^2 ) = y (^) n + h(F (y (^) n , 0) + O(h)) + τ (^) n ok
- Let u (^) n , v (^) n be the numerical solutions starting from u 0 , v 0. |u (^) n+1 − v (^) n 1 | = |u (^) n + hF (u (^) n , h) − (v (^) n + hF (v (^) n , h))| ≤ |u (^) n − v (^) n | + h|F (u (^) n , h) − F (v (^) n , h)| ≤ (1 + h L˜)|u (^) n − v (^) n |
⇒ |u (^) n − v (^) n | ≤ (1 + h L˜)n^ |u 0 − v 0 | ≤ enhL˜^ |u 0 − v 0 | = e Lt˜^ |u 0 − v 0 | ok
- define en = y (^) n − u (^) n en+1 = y (^) n+1 − u (^) n+1 = y (^) n + hF (y (^) n , h) + τ (^) n − (u (^) n + hF (u (^) n , h)) = · · ·
|en+1 | ≤ (1 + h L˜)|en | + τ (^) n ,... |en | ≤ e Lt˜^ |e 0 | + e^
Lt˜ (^) − 1 h L˜ ·^ τ^ ,^ τ^ = max^ τ^ n^ ok
note
- A method of the form u (^) n+1 = u (^) n + hF (u (^) n , h) is an explicit 1-step method. This includes Runge-Kutta and Taylor series methods.
- If τ = O(hp+1^ ), then |u (^) n − y (^) n | = O(hp^ ), i.e. the global error is one order lower than the local truncation error. We say that the method is p-th order accurate. generalization
- systems y ′ 1 =. f 1 (y 1 ,... , y (^) N ) .. y ′ N = fN (y 1 ,... , y (^) N )
⇒^ y^ ′^ =^ f^ (y)^ ,^ y^ =
^ y ...^1 y (^) N
, f =
^ f ...^1 fN
- non-autonomous equations
y ′^ = f (y, t) , yy^1 2 ==^ yt
} ⇒
( (^) y (^1) y (^2)
)′
( (^) f (y 1 , y 2 ) 1
)
- higher order equations
y ′′^ = f (y, y ′^ ) , (^) yy 2 1 == yy ′
} ⇒
( (^) y (^1) y (^2)
)′
( (^) y (^2) f (y 1 , y 2 )
)
def : A vector norm ||y|| has the following properties.
- ||y|| ≥ 0 , ||y|| = 0 ⇔ y = 0
- ||αy|| = |α| · ||y||
- ||y + u|| ≤ ||y|| + ||u|| ex ||y|| 1 = ∑N i=1^ |y^ i^ |^ ,^ ||y||∞^ = max^ |y^ i^ |^ ,^ ||y||^2 =
∑^ N i=1^ |y^ i^ |
2
1 / 2
note Given a vector norm ||y||, the subordinate matrix norm is ||A|| = max y+ =0^ ||||Ayy||||. ||A|| satisfies properties 1, 2, 3 above and 4, 5 below.
- ||Ay|| ≤ ||A|| · ||y|| , 5. ||AB|| ≤ ||A|| · ||B||
ex y ′^ = Ay : linear system assume A is diagonalizable : A = XDX −^1 D = diag(λ 1 ,... , λ (^) N ) , X = [x 1 · · · x (^) N ] , Ax (^) j = λ (^) j x (^) j , j = 1,... , N y(t) = α 1 eλ^1 t^ x 1 + · · · + αN eλ^ N^ t^ x (^) N u (^) n+1 = u (^) n + hAu (^) n = (I + hA)u (^) n u (^) n = (I + hA)n^ u 0 = α 1 (1 + hλ 1 )n^ x 1 + · · · + αN (1 + hλN )n^ x (^) N absolute stability ⇔ |1 + hλj | ≤ 1 , j = 1,... , N Hence to ensure absolute stability for a diagonalizable linear system y ′^ = Ay, it is sufficient to ensure absolute stability for the scalar test equation y ′^ = λy for each eigenvalue of A. ex ( (^) y (^1) y (^2)
)′
) ( (^) y (^1) y (^2)
) ⇒ λ 1 = − 20 , λ 2 = − 2 : stiff system
y(t) = α 1 e−^20 t^ x 1 + α 2 e−^2 t^ x 2 : fast component + slow component u (^) n = α 1 (1 − 20 h)n^ x 1 + α 1 (1 − 2 h)n^ x (^2) absolute stability ⇔ | 1 − 20 h| ≤ 1 , | 1 − 2 h| ≤ 1 ⇔ − 1 ≤ 1 − 20 h ≤ 1 , − 1 ≤ 1 − 2 h ≤ 1 ⇔ − 1 ≤ 1 − 20 h , − 1 ≤ 1 − 2 h ⇔ h ≤ 0. 1 , h ≤ 1 Hence we must choose h ≤ 0 .1 to ensure that Euler’s method is absolutely stable. note The fast component of the exact solution decays rapidly, and after some time, essentially only the slow component is present, i.e. y(t) ∼ α 2 e−^2 t^ x 2. The slow component by itself requires only h ≤ 1 for absolute stability, but we must choose h ≤ 0 .1 to ensure absolute stability of Euler’s method applied to the system.
ex :
( (^) y (^1) y (^2)
)′
) ( (^) y (^1) y (^2)
) ,
( (^) y (^1) y (^2)
) 0 =
) , solid line = y 1 (t)
! (^30) 0.5 1 1.5 2
! 2
! 1
0
1
2
3 forward Euler , h = 0.1025 , not AS
! (^30) 0.5 1 1.5 2
! 2
! 1
0
1
2
3 forward Euler , h = 0.1 , AS
! (^30) 0.5 1 1.5 2
! 2
! 1
0
1
2
3 forward Euler , h = 0.09 , AS
! (^30) 0.5 1 1.5 2
! 2
! 1
0
1
2
3 backward Euler , h = 0.1025 , AS
! (^30) 0.5 1 1.5 2
! 2
! 1
0
1
2
3 backward Euler , h = 0.1 , AS
! (^30) 0.5 1 1.5 2
! 2
! 1
0
1
2
3 backward Euler , h = 0.09 , AS
region of absolute stability for 4th order Runge-Kutta k 1 = f (u (^) n ) , k 2 = f (u (^) n + h 2 k 1 ) , k 3 = f (u (^) n + h 2 k 2 ) , k 4 = f (u (^) n + hk 3 ) u (^) n+1 = u (^) u + h 6 (k 1 + 2k 2 + 2k 3 + k 4 )
y ′^ = λy ⇒ u (^) n+1 =
1 + hλ + h^22 λ 2 + h^36 λ 3 + h 244 λ^4
(^) u (^) n
!! (^45)! 4! 3! 2! 1 0 1 2 3
! 3
! 2
! 1
0
1
2
3
4
The region of absolute stability for RK4 : { (^) 1. is larger than the region for forward Euler,
- contains an interval on the imaginary axis.
multistep methods Adams-Bashforth u (^) n+1 = u (^) n + h(β 0 f (u (^) n ) + β 1 f (u (^) n− 1 ) + · · · + βk f (u (^) n−k )) : (k + 1)-step method , explicit idea
y ′^ = f (y) ⇒ y (^) n+1 = y (^) n + ∫ (^) t (^) n+ t (^) n^ f^ (y(t))^ dt approximate f (y(t)) by an interpolating polynomial p(t) p(tn ) = f (u (^) n ) p(tn.− 1 ) = f (u (^) n− 1 ) .. p(tn−k ) = f (u (^) n−k ) def ∇fn = fn − fn− 1 : backward difference ∇ 2 fn = ∇(∇fn ) = ∇(fn − fn− 1 ) = ∇fn − ∇fn− 1 = fn − fn− 1 − (fn− 1 − fn− 2 ) = fn − 2 fn− 1 + fn− 2 ...
If (^) n = fn , S (^) − fn = fn− 1 ⇒ ∇ = I − S (^) −
p(t) = fn + (t − tn )∇ hf n + (t − tn )(t − tn− 1 )∇^
(^2) fn 2 h^2
- · · · + (t − tn ) · · · (t − tn−k+1 )∇^
k (^) fn ︸ ︷︷ ︸^ k!^ hk k terms note p(t) is a polynomial of degree k , it resembles a Taylor polynomial
p(t) = (^) i∑=0k (t − tn ) · · · (t − tn−i+1 )∇^ i (^) fn i! hi
