Download Spring 2005 HW3 Solutions for Numerical ODEs and more Assignments Mathematics in PDF only on Docsity!
Amath-Math 586/Atm S 581 Spring 2005
Homework 3 Solutions
1a) Assume A (t) is rapidly driven into a balance between its source 1 - cos( t ) and its sink -λ A. Then A ( t ) ≈ A 1 λ−1, where A 1 ( t ) = 1 - cos( t ). This dominant balance is not exact, but by
writing A in a perturbation series A = A 1 ( t )λ−1^ +Α 2 ( t )λ−2^ + ... , substituting into
dA/dt = -λ A + 1 - cos( t ),
and matching terms in powers of λ−1, we obtain Α 2 ( t ) = - dA 1 /dt = sin( t ), etc. We now
substitute this perturbation series into the second ODE
dB/dt + B = λ A = A 1 ( t ) + Α 2 ( t )λ−1^ + ... = 1 - cos( t ) + Ο( λ−1).
Solving this ODE subject to the initial condition B (0) = 0 using the integrating factor et :
[ ] ' 1 ( ) (0) 0 1 cos( ') ' ( ) t t t B t e = B + (^) ∫ e − t dt + O λ−
' (1 ) ' 1 Re 0 ' ( )
t (^) t i t = (^) ∫ ^ e − e +^ dt + O λ−
(1 ) 1 Re 1 ( 1 ) 1
i t e t^ e O i
− ^ − + − (^) +
Now, (1 ) (^1 1 ) e Re (1 )( cos 1 sin ) (cos sin ) 1 1 2 2
i t e (^) i e t (^) t iet t et t t i
=^ ^ −^ −^ +^ ^ =^ +
R −
so B t ( ) = 1 − e − t^^ − cos^ t + sin t − e − t^^ / 2 + O ( λ−^1 ) = 1 − 0.5 e − t − 0.5 cos( t + sin t (^) )+ O ( λ−^1 )
1b-c) The script hw3p1bc.m generated the desired log-log plot (below) of solution increments ε vs. ∆ t for backward Euler method and BDF2 methods. As expected, ε is O(∆ t ) for BE and O(∆ t^2 ) for BDF2. For ∆ t > , the BE method is more accurate, but the BDF2 method becomes superior for smaller ∆ t.. Both methods were coded by writing the system in the form: d φ / d t = Mφ + s ( t ,
where
A t B t
φ = , M , and
1 cos( ) ( ) 0
t t
s.
In this format, the methods are: BE: (^) ( I − M ∆ t (^) ) φ n^ +^1 = φ n^ + ∆ t s ( tn +^1 ),
BDF2: (^) ( 3 I − 2 M ∆ t (^) ) φ n^ +^1 = 4 φ n^ − φ n^ −^1 + 2 ∆ t s ( t n +^1 ).
In this problem M is tridiagonal, so one can sequentially solve for An +1^ then Bn +1, but I just used the Matlab Gaussian elimination function M.
10 −2^10 −1^100 10 −
10 −
10 −
10 −
∆ t
ε^ B−E BDF ∆ t^1 fit line ∆ t^2 fit line
2a) The problem can be written in characteristic form as du/dt = 0 ( u = constant) along characteristics dx/dt = 1+9 x. Integrating the characteristic equation from a arbitrary starting time(0, t 0 ) at which the characteristic passes through x = 0, we have
t t (^) 0 0 '/(1 9 ') log(1 9 ) / 9.
x
− = ∫ dx + x = + x
Substituting x = 1, and setting τ = log(10)/9 = 0.256 we conclude that
u (1, t 0 ) = u (0, t 0 + τ) =
t t
^ <
Only one BC is required per characteristic, and this has been specified at x = 0, so no additional BC is required at x = 1.
2b) The stencil for the leapfrog centered-in-space method (right) includes the analytical domain of dependence (the characteristic dt/dx = c-1 through (xj,tn+1)) within the numerical domain of dependence if the Courant number c ∆ t/ ∆ x < 1. The maximum c across the domain is cmax
= ( ) , so the CFL stability limit is ∆ t ≤ ∆ x/ c
0 1
max 1 9 10 x
x ≤ ≤
2c) A second-order accurate one-sided approximation to ux at the maximum gridpoint j = N is
1 2 (1, ) 2 2 2 1.5 2 0.5 / B n B n n n n
u x tn δ x uN δ x u N u N uN uN
≈ − = − −^ + − ∆ x
This can be seen by Taylor expanding one-sided differences of span ∆ x and 2∆ x and finding a linear combination of these that cancels the O(∆ x ) part of their truncation errors:
δ xB u x ( ) = ( u x ( ) − u x ( − ∆ x ) / ∆ x = ux − ∆ x u ( xx / 2 ) + O ( ∆ x^2 )
δ 2 B x u x ( ) = ( u x ( ) − u x ( − 2 ∆ x ) /(2 ∆ x ) = ux − ∆ x u ( xx ) + O ( ∆ x^2 )
so 2 δ 2 B x^ u x ( ) − δ 2 Bx u x ( ) = ux + O ( ∆ x^2 ).
2d) The Matlab script hw3p2d.m implements the leapfrog-centered in space method with a midpoint RK2 starting step, the one-sided difference at the right boundary, and optional Asselin filtering, using ∆ x/ cmax = 0.05 and ∆ t = 0.005, producing the plots below. Asselin filtering is required to control the 2∆ t errors that otherwise dominate u (1, t ). With Asselin filtering there is a slight overshoot (dispersion error) in the solution, which is what we would expect from the even-order spatial differencing, but the main error is diffusion of the exact step solution due to the temporal smoothing of high frequencies from the Asselin time- filtering scheme.
−10 0 0.1 0.2 0.3 0.4 0.
−
0
5
10
t
u(1,t)
No Asselin filtering
−0.2 0 0.1 0.2 0.3 0.4 0. 5
0
1
t
u(1,t)
Asselin, γ = 0. Leapfrog Exact
Leapfrog Exact
