Exercise Numerical Analysis of Differential Equations

Exercise

Numerical Analysis of Differential Equations

Fakult¨ at f¨ ur Mathematik und Informatik TU Bergakademie Freiberg

Dr. M. Helm / Dipl.-Math. St. Pacholak

August 13, 2015

Contents

1 Initial Value Problems - Repetition 3

1.1 Comprehension Questions . . . 3

1.2 Mixing Tank. . . 4

1.3 Interconnected Tanks . . . 5

1.4 Even more Tanks . . . 6

1.5 Higher Order IVP . . . 7

2 Initial Value Problems - Single-step Methods 8 2.1 Comprehension Questions . . . 8

2.2 EULER Methods . . . 9

2.3 BUTCHER Arrays . . . 10

2.4 RUNGE-KUTTA Method . . . 11

2.5 Consistency . . . 12

2.6 Self study . . . 13

3 Initial Value Problems - Multi-step Methods 14 3.1 Comprehension Questions . . . 14

3.2 Characteristic Polynoms . . . 15

3.3 Order of Concistency . . . 16

3.4 Concistency, Zero-Stability & Convergence . . . 17

4 Initial Value Problems - MATLAB Training 18 4.1 Figures und Functions . . . 18

4.2 Subplots and Function Handles . . . 19

5 Boundary Value Problems - Elliptic BVP 20 5.1 Comprehension Questions . . . 20

5.2 DIRICHLET Boundaries . . . 21

5.3 Von NEUMANN Boundaries . . . 22

5.4 Universal Grid Size . . . 23

5.5 DIRICHLET & von NEUMANN Boundaries . . . 24

5.6 Self Study . . . 25

1

6 Boundary Value Problems - 1D & 2D 26 6.1 Square Grid Cells . . . 26 6.2 Rectangular Grid Cells . . . 27 6.3 Discretization in Space and Time . . . 28

7 Assignment - Restricted 3-Body Problem 29

1 Initial Value Problems - Repetition

1.1 Comprehension Questions

I) The initial value problem y⁰(x) = f(x, y(x)) with y(a) = ya has a unique solution on [a,b], if:

A f is continuous.

B f is continuous and Lipschitz continous in its second argument.

C f is differentiable.

D f is Lipschitz continous in its first argument.

II) Which of the following are NOT first order initial value problems?

A ^dy_dx = 2(x+y) ,y(0) =y₀

B y⁰⁰−2y⁰+y=te^t−t , y(0) = 0 C y⁰ = ¹₂y² , y(1) = 1

D y⁰₁(t) y⁰₂(t)

!

= 3 0

0 −2

! y₁(t) y2(t)

!

, y 0 0

!

= 4

7

!

III) Which of the following statements are correct?

A y⁰ = 3 0

0 −2

!

has the Eigenvalues λ₁ = 3 andλ₂ =−2.

B The initial value problem ^dQ_dt = ^r₄ − ₁₀₀^r Q with Q(0) = Q₀ has no solution in [0, t_e].

C ^dQ_dt = ^r₄ −₁₀₀^r Q describes a predator-prey-interaction.

3

1.2 Mixing Tank

At time t = 0 a tank contains Q₀ kg of salt dissolved in 100 m³ of water. Assume that water containing 0.25 kg/m³ of salt is entering and leaving the tank at a rate ofrm³/min.

a) What’s the according differential equation for this initial value problem (IVP)?

b) Show that the IVP has a unique solution using the Picard-Lindel¨of theorem!

c) Determine the general and exact solution of the IVP!

d) Take MATLAB or your pocket calculator to draw some pictures of the solution for several initial values Q₀ and r= 3.

1.3 Interconnected Tanks

Remembering the Mixing Tank, we extend our experiment by adding another tank and two connecting pipes between them, as shown in the figure.

a) Create the according system of differential equation for this initial value problem (IVP)?

b) Solve the IVP with the given values using eigenvalues and their according vectors to diagonalize the matrix A!

5

1.4 Even more Tanks

Considering the shown experimental setup as an extension of the Interconnected Tanks, we added flow inlets for each tank and an outlet to get a non discharging system.

a) Determine the exact solution of the linear problem with the initial valuesQ₁(0) = 55 and Q₂(0) = 26!

b) According to the laboratories resources we can’t use these tank sizes and have to work with two tanks with only 10 m³but same concentrations of salt in it. Rearrange and solve your IVP using half filled tanks as starting point!

c) Create the equation system for determining the concentration of salt in both tanks in setup b) after 30 min, when the outlet plugged up after 20 min. What do you recognize when trying to solve it?

1.5 Higher Order IVP

Rewrite the following higher order initial value problems from science and technical application as first order systems!

a)y⁰⁰−2y⁰+y=te^t−t 0≤t ≤1 y(0) = 0, y⁰(0) = 0

b)y⁰⁰⁰ = ¹₂y² t≥1 y(1) = 1, y⁰(1) =−1, y⁰⁰(1) = 2 c)t²y⁰⁰(t)−2ty⁰(t) + 2y(t) =t³lnt t≥1 y(1) = 1, y⁰(1) = 0

7

2 Initial Value Problems - Single-step Methods

2.1 Comprehension Questions

I) Which of the following BUTCHER arrays describe an explicit method?

A) 0 B) 1 1

1/2 1/2 1/2 1/4 1/4

1/2 1/4 1/4 1 1/4 1/4 1/2

1 0 -1 2 1 0 -2 2 1

1/6 0 2/3 1/6 1/6 0 2/3 1/6

C) 0 -1 2 0 D) A

1/2 1/4 1/4 B

1/2 1/2 1/3 5/12 -1/12 C

1 1 3/4 1/4 D

1/6 0 2/3 1/6 3/4 1/4

II) Which of the following methods are consistent?

A) 0 B) 0

1/2 1/2 1/2 1/2

1/2 1/6 1/2 1/2 -1/4 1/4

1 0 -1 2 1 0 -1 2

1/6 0 2/3 1/6 1/6 0 2/3 1/6

C) 0

1/2 1/2 A

1/2 1/4 1/4 B

1 0 -1 2 C

1/6 0 2/3 1/6

III) Which method is described through the Butcher array 0

1 ? A Euler modified

B Euler explicit

C Euler improved

D Euler implicit

2.2 EULER Methods

Determine an approximated solution for the initial value problem y⁰ =−(y+ 1)(y+ 3), 0≤t ≤2 y(0) = -2

applying

a) the explicit EULER scheme, b) the modified EULER method, c) the improved EULER method

with step sizeh= 0.5 andh= 0.2. Compare your approximations with the corresponding values of the exact solution y(t) =−3 + 2(1 +e^−2t)⁻¹.

Hint: For the improved EULER method take y_n+1 = y_n + ^h₂[f(t_n, y_n) + f(t_n +h, y_n + hf(t_n, yn))]

9

2.3 BUTCHER Arrays

Rewrite each of the following BUTCHER/RUNGE-KUTTA arrays as a set of explicit formulae for the calculation of y_j+1. Do they describe explicit or implicit methods?

a) four-stage ENGLAND formula 0

1/2 1/2 1/2 1/4 1/4

1 0 -1 2

1/6 0 2/3 1/6 b) BUTCHER’s formula

0

1/8 1/8

1/4 0 1/4

1/2 1/2 -1 1

3/4 3/16 0 0 9/16

1 -5/7 4/7 12/7 -12/7 8/7

7/90 0 32/90 12/90 32/90 7/90 c) (one) RADAU IIA method

1/3 5/12 -1/12

1 3/4 1/4

3/4 1/4

2.4 RUNGE-KUTTA Method

For the solution y(x) of the initial value problem

y⁰ =−(y+ 1)(y+ 3), 0≤t ≤2 y(0) = -2

find an approximation of y(1) using the classical RUNGE-KUTTA method with step size h = 0.5. Compare with the solutions in 2.2.

11

2.5 Consistency

Check the concistency order of HEUN’s method by verification of the order conditions (for explicit RUNGE-KUTTA methods up to order three) known from the lecture.

0

1/3 1/3

2/3 0 2/3 0

1/4 0 3/4

2.6 Self study

Determine an approximated solution to the following system of ODEs.

˙y(t) = 6 −3 2 1

!

y(t), y(0) = 5 3

!

, 0≤t≤1.

Therefore apply

a) the explicit EULER method with step sizeh = 0.25, b) HEUN’s method with step sizeh = 0.5.

Compare the approximations with the exact solution y(t) = exp(3t) 3exp(4t)

exp(3t) 2exp(4t)

! −1 2

!

Hint: You have to use vector versions of both methods, for instance y_n+1 = y_n+hf(t_n,y_n) and y₀ =y(0) in the EULER case.

13

3 Initial Value Problems - Multi-step Methods

3.1 Comprehension Questions

I) Which of the following statements are correct? A multi-step method is:

A convergent, if it’s consistent and zero-stable.

B consistent, if it’s convergent and zero-stable.

C zero-stable, if it’s convergent and consistent.

II) Determine the non-consistent methods!

A y_n+1 =y_n

B y_n+3 =−y_n+2+y_n+1+y_n+ 2h(f(y_n+2) +f(y_n+1)) C y_n+1−y_n−1 = 2hf(t_n, y_n)

D y_n+2 =y_n+^h₃(f_n+ 4f_n+1+f_n+2)

III) Which of the following methods are no real multi-step methods?

A y_n+1 =y_n

B y_n+3 =−y_n+2+y_n+1+y_n+ 2h(f(y_n+2) +f(y_n+1)) C y_n+1−yn−1 = 2hf(t_n, y_n)

D y_n+2 =y_n+^h₃(f_n+ 4f_n+1+f_n+2)

3.2 Characteristic Polynoms

Determine the characteristic polynomialsρ(ζ) andσ(ζ) for the following linear multistep methods. Verify the conditions

Pr j=0

αj = 0 und ^Pr

j=0

jαj = ^Pr

j=0

βj

for

a) the 3-step ADAMS-BASHFORTH method, b) the 3-step ADAMS-MOULTON method, c) 2-step SIMPSON method.

15

3.3 Order of Concistency

Determine the order of consistency for

a) the explicit NYSTR ¨OM method withr = 2 (midpoint rule) yn+1−yn−1 = 2hf(t_n, yn),

b) the implicit NYSTR ¨OM method with r = 2 (SIMPSON rule)

y_n+2−y_n= ^h₃{f(t_n, y_n) + 4f(t_n+1, y_n+1) +f(t_n+2, y_n+2)}.

3.4 Concistency, Zero-Stability & Convergence

Which of the following linear multistep methods are convergent? If a method is not con- vergent: is it consistent, not zero-stable, or both?

a) yⁿ⁺² = ¹₂yⁿ⁺¹+ ¹₂yⁿ+ 2hf(yⁿ⁺¹) b) yⁿ⁺¹ = yⁿ

c) yⁿ⁺⁴ = yⁿ+ ⁴₃h{f(yⁿ⁺³) +f(yⁿ⁺²) +f(yⁿ⁺¹)}

d) yⁿ⁺³ =−yⁿ⁺²+ yⁿ⁺¹+ yⁿ+ 2h{f(yⁿ⁺²) +f(yⁿ⁺¹)}.

17

4 Initial Value Problems - MATLAB Training

Hint: This exercises should be done during the problem session in the PC-Lab MIB 2.10.

Label the axes, titles and legend of the subplots comprehensibly.

4.1 Figures und Functions

Consider the initial value problem

y⁰(t) = tsin(2y), y(0) = ^π₄. a) Approximate the exact solution

y(t) = arctan(e^t2)

on the time interval [0,3] with the modified EULER and HEUN‘s third-order method using self-defined functions. Draw figures of the exact and the numerical solutions for three different step sizes h<1.

b) Approximate the solutions with the explicit fourth-order ADAMS-BASHFORTH method. Take the first values of HEUN‘s method in a) as a startup and add the curve into the 3 figures.

c) Draw a figure with logarithmic axes of the error maxi|yi −y(ti)| of the numerical solutions for step sizes between 1 and 10⁻⁴.

4.2 Subplots and Function Handles

The ODE

y⁰(t) =λ(y(t)−g(t)) +g⁰(t) has the general solution

y(t) = Ce^λt+g(t), C ∈R

which, for Re λ < 0, consists on a decaying transientCe^λt and a steady state component g(t). We consider the case g(t) = arctan t, λ = -10.

a) Draw a subplot of the general solution, i.e. plot the graph ofy(t) for several choices of C ∈[-5,5].

b) Find the exact solution for the initial conditiony(0) = 0.

c) Try to approximate the solution of this IVP with the explicit EULER method using function handles. Plot the numerical solutions on the time interval [0,5] forh= 0.5, h = 0.25 and h = 0.1 in different subplots.

d) Repeat the experiment from (c) with the implicit EULER method using function handles and symbolic functions. Comment your observation.

19

5 Boundary Value Problems - Elliptic BVP

5.1 Comprehension Questions

I) Which of the following matrices describes the BVP -u⁰⁰(x) = 2, x ∈ Ω := (0,1), u(0) = 0, u⁰(1) = 0?

A 2 −1

−1 2

! u1

u₂

!

= −2− ₂₄₃⁵

−1− ₂₄₃⁴⁰

!

B







2 −1 0

−1 2 −1

1 −1 1













u₁ u₂ u₃





 =







−2− ₂₄₃⁵

−1− ₂₄₃⁴⁰

1 3







C







2 −1 0

−1 2 −1

0 −1 1













u1 u2 u3





 =







−2− ₂₄₃⁵

−1− ₂₄₃⁴⁰

2 3







II) Which of the following statements are wrong for the BVP -u⁰⁰(x) = 2, x ∈ Ω := (0,1), u(0) = 0, u⁰(1) = 0?

A The BVP has two DIRICHLET boundary conditions!

B The BVP has a von NEUMANN boundary condition!

C The BVP is two-dimensional!

D The BVP is one-dimensional!

III) When discretizing a two-dimensional BVP the 5-point stencil is used.

Which statement is NOT correct?

A It’s also called 5-point discrete Laplacian!

B It can be derived from TAYLOR expansion forumula!

C It contains the 6 neighbors of each discretized point!

D This equation is not valid for boundary points!

5.2 DIRICHLET Boundaries

Consider the following elliptic boundary value problem of second order in the one dimen- sional space

-u⁰⁰(x) = -5x³ x ∈ Ω = (0,1) u(0) = -2;

u(1) = -1

a) Check that u(x) = ¹₄x⁵+³₄x - 2 is the solution to the above boundary value problem.

b) Discretize the boundary value problem with central finite differences on a uniform grid with mesh width ∆x = h = ¹₃.

Write down the associated linear system and determine a numerical approximation to the solution of the boundary value problem by solving this linear system.

c) Compare your numerical solution with the exact one.

21

5.3 Von NEUMANN Boundaries

Modify the boundary value problem from Exercise 5.2 to -u⁰⁰(x) = -5x³ x ∈ Ω = (0,1) u(0) = -2;

u⁰(1) = 2

a) Show, that u(x) = ¹₄ x⁵ + ³₄ x - 2 is still a solution of this boundary value problem.

b) Discretize the boundary value problem with central finite differences on a uniform grid with mesh width ∆x = h = ¹₃.

Write down the associated linear system and determine a numerical approximation to the solution of the boundary value problem by solving this linear system.

c) Compare your numerical solution with the exact one.

5.4 Universal Grid Size

Generate the matrices and right hand sides of Ax = b, when an arbitrary number of grid points is used. Try to simplify these matrices through (n+1) block structures?

23

5.5 DIRICHLET & von NEUMANN Boundaries

Consider the following second order elliptic boundary value problem in the one dimen- sional space

-u⁰⁰(x) = 2; x∈Ω := (0,1), u(0) = 0;

u⁰(1) = 0;

a) Verify that the solution of the boundary value problem is given byu(x) =x(2−x).

b) Discretize the boundary value problem using central differences on a uniform grid with mesh width ∆x = ¹₃. Write down the associated linear system and determine a numerical approximation to the solution of the boundary value problem by solving this linear system.

c) Compare your numerical solution with the exact one. What happens if the NEU- MANN boundary condition is replaced by u(1) = 0? Explain your observations.

5.6 Self Study

Determine the approximate stationary temperatur distribution in a thin quadratic metal plate with a side length of 0.5 m. Two adjacent are hold on a temperature of 0^◦C. On the other boundaries are hold on a temperature of 0^◦C to 100^◦C.

(Hint: Ifxandyare the space coordinates, the problem can be describe by the stationary heat equation ∆u(x,y) = 0 together with appropriate boundary value conditions. for the solution of the linear system you should use software like MATLAB. for purposes of comparison notice, that the exact solution is given by u(x, y) = 400xy if the boundaries with zero boundary condition are placed along the coordinate axes.)

25

6 Boundary Value Problems - 1D & 2D

6.1 Square Grid Cells

Discretize the stationary heat equation

-∆u(x, y) = -3 ((x,y) ∈Ω) on the following mesh:

a) Use u = 0 as DIRICHLET boundary values on the open circles and u = 1 on the open squares.

b) Replace the Dirichlet boundary condition on the left and on the bottom-right by

δu

δn = 0 on the lines 1 3 and 17 19.

6.2 Rectangular Grid Cells

A rectangular silver plate of 6 cm x 5 cm is heated uniformly at each point with a constant energy rateq= 6.2802 W/cm³. Assume that the temperature along the boundaries is hold at

u(x,0) =x(6−x), u(x,5) = 0 (0< x <6), u(0, y) = y(5−y), u(6, y) = 0 (0< y <5).

Then, the stationary temperature distribution can be described by the POISSON equation

∆u(x, y) =−_K^q for (x, y)∈Ω := (0,6) × (0,5),

together with the above boundary conditions. The constant K = 4.3543 W/(cm· K) is the thermal conductivity.

Approximate the temperature distribution by using the finite difference method with mesh widths ∆x = 0.4 and ∆y = ¹₃. Try also some other mesh widths.

27

6.3 Discretization in Space and Time

Consider the initial-boundary value problem

u_t(x, y)− ¹₂u_xx(x, t) = 0 in Ω = (0,1) x (0,T), u(x,=) =x(1−x) f¨urx∈ (0,1),

u(0, t) =u(1, t) = 0 f¨urt ∈(0,T) for an unknown function u = u(x;t).

a) Approximate u= u(x;t) at time t = 0.1 with an explicit difference scheme. Take a step size of ∆x = 0.25 in space and ∆t = 0.1 in time.

b) Approximate u = u(x, t) at time t = 0.1 with a purley implicit difference scheme.

Take the same step size as in a).

c) For which step sizes ∆t in time is the explicit and the purely implicit difference scheme numerically stable, if a step size ∆x = 0.25 in space is provided?

7 Assignment - Restricted 3-Body Problem

We consider a satellite in the gravity field of earth and moon. We assume that the move- ment of the three bodies takes place within a fixed plane, and that the two heavy bodies are rotating around their barycenter with constant distance and with constant angular velocity. This means, the influence of the satellite on the orbits of earth and moon can be omitted - that’s why this is called a restricted 3-body problem.

With respect to a coordinate system which is coupled to this rotation (i.e. earth and moon have fixed positions in this system), the satellite’s orbit (x, y) = (x(t), y(t)) can be described by a system of second order differential equations:

x⁰⁰ =x+ 2y⁰−µ_[(x+y)^x+µ2+y²]^3/2 −µ_[(x−µ)^x−µ2+y²]^3/2

y⁰⁰ =y−2x⁰−µ_[(x+µ)2^y_+y2]3/2 −µ_[(x−µ)2^y_+y2]3/2

Hereµ= 1/82.45 is the ration between the mass of the moon and the cumulative mass of earth and moon, andµ= 1 - µ. The length unit is the distance between earth and moon.

The earth is placed in the origin and the moon on the positive real axis. The time unit is chosen in a way, that the angular velocity of the rotation is one, i.e. the moon needs 2π time units for one orbit.

a) Rewrite the second order differential equation system as a first order system. Apply the substitutions y₁(t) =x(t), y₂(t) = y(t) and y₃(t) =x⁰(t), y₄(t) =y⁰(t).

b) Solve the first order ODE system using the MATLAB routine [t,y] = ode45(@f,y0,[0 tend],options)

with the following initial values y0 in t = 0:

y1(0) = 1.2, y₂(0) = 0, y₃(0) = 0, y₄(0) = -1.04935750983031990726 Therefore chosen and appropriate end tend for the integration interval.

29

c) Specify or control via the fourth input parameter options:

- the tolerances RelTol and AbsTol for the stepsize control ofode45, - the plot function for the phase plot (y₁, y₂) (corresponds to (x(t), y(t))), - the detection of local minima and maxima of the orbit/trajectory (points with

horizontal or vertical tangent) by use of an event function.

Hint: Analyze the additional output parameters in the statement [t,y,te,ye,ie] = ode45(@f,y0,[0 tend],options).

d) Observe what happens if

- tend will be enlarged or reduced,

- the fourth component in the initial condition y₄(0) is rounded (for instance) to -1.05

e) Modify the event function so that ode45 terminates after the calculation of a closed trajectory.

f) Solve the 3-body problem with a classical RUNGE-KUTTA method of fourth order an constant stepsize h= (t_end−t₀)/n. How many integration steps n are needed to get a similar trajectory as withode45?

Use the plot statement to generate a picture of the orbital trajectory.

g) Solve the 3-body problem a second time with ode45 using the following data µ= 0.012277471,

y₁(0) = 0.994, y₂(0) = 0, y₃(0) = 0, y₄(0) = -2.00158510637908252240537862224.

Determine a value for tend resulting in a closed trajectory. How many integration steps are needed for the associated integration interval?

Hint: The number of integration steps in ode45 depends on the chosen tolerances RelTol and AbsTol.

h) Try to reproduce the second orbit with the classical RUNGE-KUTTA method with fixed stepsize.