Finite-Difference Solver - Numerical Methods for Partial Differential Equations - Lecture Slides, Slides of Mathematical Methods for Numerical Analysis and Optimization

Download Finite-Difference Solver - Numerical Methods for Partial Differential Equations - Lecture Slides and more Slides Mathematical Methods for Numerical Analysis and Optimization in PDF only on Docsity!

Numerical Methods for Partial

Differential Equations

Homework 3

Q1) Build a finite-difference solver for Q1a) use your Cash-Karp Runge-Kutta time integrator from HW2 fortime stepping

Q1b) use the 4th^ order central difference in space (periodic domain) Q1c) perform a stability analysis for the time-stepping (based on visual inspection of the CKR-K stability region containing the imaginary axis)

Q1d) bound the spectral radius of the spatial operator Q1e) choose a dt well in the stability region Q1f) perform four runs with initial condition (use M=20,40,80,160) and compute maximum error at t=

Q1g) estimate the accuracy order of the solution. Q1h) extra credit: perform adaptive time-stepping to keep the local truncation error fromtime stepping bounded by a tolerance.

u c u

t x

( ) (^ )^ [ ] 4cos^2

u x ,0 = e −^ π^ x , x ∈ 0,

Where the b,c coefficients are

http://www.library.cornell.edu/nr/bookcpdf/c16-2.pdf

c= 37/378 [ 0 250/621 125/594 0 512/1771]

b

= ^ 

Homework 3

Q1) Build a finite-difference solver for Q1a) use your Cash-Karp Runge-Kutta time integrator from HW2 for time stepping Q1b) use the 4th^ order central difference in space (periodic domain) Q1c) perform a stability analysis for the time-stepping (based on visual inspection of the CKR-K stability region containing the imaginary axis)

Q1d) bound the spectral radius of the spatial operator Q1e) choose a dt well in the stability region Q1f) perform four runs with initial condition (use M=20,40,80,160) and compute maximum error at t=

Q1g) estimate the accuracy order of the solution. Q1h) extra credit: perform adaptive time-stepping to keep the local truncation error fromtime stepping bounded by a tolerance.

u c u

t x

( ) (^ )^ [ ] 4cos^2

u x ,0 = e −^ π^ x , x ∈ 0,

Code

• Each time step consists of evaluating the 6 intermediate vectors k1,k2,..,k
• And piecing them together.
• Notice – here I set up u at the start as a vector, and delta4 accepts a vector (and dx) as arguments and returns a vector.
• i.e. each of k1,k2,..,k6 is a vector.

Alternate Code

• I could also have used loops

• Note: u is a row vector of length M so the k’s are

a matrix of size 6 x M

Cash-Karp

• The Cash-Karp Runge-Kutta integrator in our case is:
• Notice, we do not need the time component on the right hand side, because our finite-difference right hand side is independent of time.

1 2 21 1 3 31 1 32 2 4 41 1 42 2 43 3 5 51 1 52 2 53 3 54 4 6 61 1 62 2 63 3 64 4 65 5 1 1 1 2 2 3 3 4 4 5 5 6 6

n n n n n n

n n

k dtf u k dtf u b k k dtf u b k b k k dtf u b k b k b k k dtf u b k b k b k b k k dtf u b k b k b k b k b k u + u c k c k c k c k c k c k

Linear Stability Analysis

• We achieve the linear stability analysis by

assuming f is linear in u:

f ( u )= μ u

1 2 21 1 3 31 1 32 2 4 41 1 42 2 43 3 5 51 1 52 2 53 3 54 4 6 61 1 62 2 63 3 64 4 65 5 1 1 1 2 2 3 3 4 4 5 5 6 6

n n n n n n

n n

k dtu k dt u b k k dt u b k b k k dt u b k b k b k k dt u b k b k b k b k k dt u b k b k b k b k b k u u c k c k c k c k c k c k

cont

• We wish to remove all the intermediate variables

and express the scheme in a one-step form like:

• Where the multiplier function pi(nu) is a 6th^ order

polynomial in nu.

• It is certainly possible to find all the coefficients

of this polynomial by hand but a little messy.

• We can take a short cut – assuming Cash-Karp is

actually a 5th^ order scheme as claimed we know

that the first 6 terms of the polynomial multiplier

must be the first 6 terms of the Taylor series for

u n^ +^1 = π ν ( ) un

e^ ν

cont

• So we know that:

• Where C is to be determined.

• However, looking at the definition of the stages

we see that there is only one contribution to the

6 th^ order term:

• So

π ν ( ) = 1 + 1!^1 ν 1 + 2!^1 ν 2 + 3!^1 ν 3 + 4!^1 ν 4 + 5!^1 ν 5 + C ν^6

c b b b b b 6 65 54 43 32 21 ν^6

π ν ( ) = 1 + 1!^1 ν 1 + 2!^1 ν 2 + 3!^1 ν 3 + 4!^1 ν 4 + 5!^1 ν 5 + c b b b b b 6 65 54 43 32 21 ν^6

Matlab roots Command

• The matlab roots command can be used to find

roots of a polynomial.

• If the polynomial is say:

• Then construct a vector of coefficients (highest

order first):

• And invoke: roots ( a )

• A vector of the roots (if any) of the polynomial

will be returned.

π ( x^ ) = a 1^ + a x 2^ + a x 3 2

a = [ a 3 , a 2 , a 1 ]

Adapting RKabstab.m

• So I took the routine from the class web page and modified it slightly
• I hard coded it to use the multiplier polynomial for the Cash-Karp
• I changed the plotting to use a scatter plot

On The Imaginary Axis

Here I used Matlab’smultiplier polynomial for nu on the imaginary axis. polyval to evaluate the Conclusion – the absolute stability region is justto the left of the imaginary axis (not a big issue here)

Homework 3

Q1) Build a finite-difference solver for Q1a) use your Cash-Karp Runge-Kutta time integrator from HW2 for time stepping Q1b) use the 4th^ order central difference in space (periodic domain) Q1c) perform a stability analysis for the time-stepping (based on visual inspection of the CKR-K stability region containing the imaginary axis)

Q1d) bound the spectral radius of the spatial operator Q1e) choose a dt well in the stability region Q1f) perform four runs with initial condition (use M=20,40,80,160) and compute maximum error at t=

Q1g) estimate the accuracy order of the solution. Q1h) extra credit: perform adaptive time-stepping to keep the local truncation error fromtime stepping bounded by a tolerance.

u c u

t x

( ) (^ )^ [ ] 4cos^2

u x ,0 = e −^ π^ x , x ∈ 0,