Regulations:

● Python is the only computer language to be used in this course

● You are not allowed to use high level Python functions (if not explicitly mentioned), such as diff,  int,  taylor,  polyfit, or ODE solvers or corresponding functions from scipy etc.  You can of course use mathematical functions such as exp, ln, and sin and make use of numpy .array.

● In your report, for each question you should add the Python code, and worked-out results, including figures or tables where appropriate.  Use either a Python script and a pdf for the report or combine everything in a Jupyter notebook.  We need to be able to run your code.

● Output numbers in floating point notation (using for example the npPrint function from the quiz).

● A concisely written report, compiled with clear presentation and well structured coding, will gain marks.

● Do not put your name but put your student identity number on the report.

● The first part of each assignment needs to be done individually by each student.

● For the second part submission in pairs is allowed, in which case one student should submit the second part and both need to indicate at the top of the their pdf for the first part the student id of their partner. Both will then receive the same mark for the second part of the assignment.

PART 1

Q 1.1.  [3/20] Consider the following model for a fish population

B\ (t) = rB(t) ╱ 1 - 、 - β

where r, K, β, α are given constants. The model assumes growth according to the logistic equation where K is the carrying capacity and the second term describes a loss rate due to predators.  It approaches a constant loss rate β for a large fish population (B - o). For a small fish population the loss is quadratic in the size of the population.                   We assume that the dimensions for B and t are fish and week, respectively. So [B] = fish and [t] = week.

(1) What dimensions do the constants r, K, β, α have to have for the model satisfies dimensional homogeneity?

(2) Define non-dimensional τ = t/T and x(τ ) = B(Tτ )/F and derive the model

x2      

x\   = x(1 - x) - Π 1

by chosing suitable characteristic time and fish scales T, F .  Check that the new parameters Π 1 , Π2  are dimensionless.

(3) Using a different choice for T, F derive the model

x2     

1 + x2   ,

with non-dimensionless (check) parameters Π3 , Π4 .

Q 1.2.  [2/20] Consider the homogeneous ODE y\  = f(y) with ∂yf(y) invertible for all y . The stiff Taylor scheme is an explicit and surprisingly enough L-stable method of the form

(*)                       yn+1 = yn + ╱ exp(h∂yf(yn)) - Im∶ ╱∂yf(yn)、-1 f(yn) .

Here we use the matrix exponential exp(A) and Im  e Rmam  is the identity matrix.      Remark: we will look at a Runge-Kutta type version of this idea in the second part.

(1) For a given tn  we can rewrite our ODE in the form

y\ (t) = F (y(t)) := Dny(t) + gn(y(t)) , y(0) = Y (tn)

with Dn  := ∂yf(Y (tn)) and gn(y) := f(y) - Dny .  Note that F (y) = f(y) so this problem still has the same exact solution Y which can be written in the form

t

Y (t) = exp(Dn(t - tn))Y (tn) +      exp(Dn(t - τ ))gn(Y (τ )) dτ .

tn

Derive the method (*) by using the approximation gn(Y (τ )) ≈ gn(Y (tn)) and as usual replacing Y (tn), Y (tn+1) by yn, yn+1, respectively.

(2) Show that the method is exact for any linear system of ODE y\  = Ay  (so the method is L-stable).

Q 1.3.  [5/20] Consider an ODE y\ (t) = f(y(t)).  An invariant (or first integral ) of this ODE is a function H : Rm  - R that satisfies VH(y) · f(y) = 0 for all y e Rm . From the definition it follows that if Y is the solution to the ODE then H(Y (t)) = 0 using chain rule and therefore H(Y (t)) = H(Y (0)) for all t (compare with the discussion of conserved quantities and invariances for kinetic equations and the Hamiltonian systems).                 Consider in the following a s-stage Runge-Kutta method with consistency order of p 2 1 given by a Butcher tableau consisting of α, β, γ .

(1) Show that the Runge-Kutta method conserves any linear invariant, i.e., for in- variants of the form H(y) = c · y with some fixed c e Rm  show that for all n: H(yn) = H(y0 ).

(2) Prove that under the condition

γiβiaj + γjβjai = γiγj

for i, j = 1, . . . , s the Runge-Kutta method maintains quadratic invariances of the form H(y) = y · Cy with a matrix C e Rmam .

Hint:

Consider H(z) = z · Cz which implies VH(z) = (C +CT )z so that H is an invariant if (C + CT)z · f (z) = z · (C + CT)f (z) = 0 for any vector z .

Use that for any vector a, b we have (a +b) · C(a +b) = a · Ca +a · (C +CT)b +b · Cb to show that the remainder R in H(yn+1) = H(yn) + R is

R = h      γiyn · (C + CT)ki + h2         γiγjkj  · Cki  .

i                                                          ij

Now you need to prove R = 0 under the given condition on β, γ which follows for example from using that Ynai  · (C + CT)ki  = Ynai  · (C + CT)f (Ynai) = 0 with Ynai = yn+h    j(s)=1 βijkj which allows us to replace yn · (C +CT)ki in the expression for R.

(3) Use the above to find all values for  d so that the Runge-Kutta method given d        d           0

by the Butcher tableau   1 - d  1 - 2d   d  maintains quadratic invariances, i.e.,

   1                  1

invariances of the form H(y) = y · Cy .

PART 2

Hand in one program listing which covers all the questions given below in a single scrip- t/notebook reusing as much as possible. Avoid any code duplications (or explain why you decide to duplicate some part of the code).

Important: In your brief discussion of your result (a paragraph with a plot or a table is sufficient) refer back to the theoretical results discussed in the lecture or from assignments.

Q 2.0.  [see online quiz on moodle]

Implement a function to performs one step of a general diagonally implicit Runge-Kutta method.   The two vectors and matrix defining the DIRK should be passed in to the function and the function needs to work for vector valued ODEs.

Once you have implemented a function of the form def  y  =  dirk(  f,Df,  t0,y0,  h,    alpha,beta,gamma) you can define a stepper  fixing alpha,beta,gamma which you can    then use in your evolve method, e.g., given the vectors and matrix alphaCN,betaCN,gammaCN defining the Crank-Nicholson method  (provided in the lecture handout) the stepper is    CrankNicholson  =  lambda  f,Df,t0,y0:    dirk(f,Df,t0,y0,alphaCN,betaCN,gammaCN) The rest of your code should then work without change. You can test your code by compar-

ing the results you get with your old steppers (forwardEuler, backwardEuler,Q11,crankNicholson). The errors should be (close to) identical.

Q 2.1.  [2/20]

Test your implementation of the Diagonally Implicit Runge-Kutta methods using our test ODE from the previous assignments: y\ (t) = F (y(t)) and y(0) = y0  with

f (t, y1 , y2 ) = 、y2 (1 y2- 2y1 ) _

for t e [0, 10] and y0 = (2, -2)T .

Compute the maximum errors and EOCs using the time step sequence  hi   =    for i = 0, . . . , 9 using N0 = 25 for the following three DIRK methods:

d        d           0

2                 2

First check the order conditions for all three method to determin if the methods have consistency order equal to 1, equal to 2, or 2 3.  Then discuss if the numerical results match your expectations.

Q 2.2.  [3/20]

Note:  Q2.3 is an extension of Q2.2 and you should combine both into one report if you can - Q2.3 includes the experiments done here so feel free to skip this question.  I just added this question for anyone who wants to only do the easier part of Q2.3.

Consider the pendulum modelled by the Hamiltonian H(x, p) =  p2  - ω0(2) cos(x).  The problem has the exact solution:

We set ω0 = 4 and simulate on the interval [0, 100].

Compute errors and EOCs for the two implicit two stage DIRK methods from the previous question using the sequence Ni  = 500 2i  for the number of time steps with i = 0, . . . , 7. In addition to the usual error vs.  h plots and discussing the EOCs also discuss how well each method conserves the Hamiltonian at least for some time steps sizes, i.e., for i = 1 and i = 3.

Carry out experiments with initial conditions (x0 , p0 ) = (0.05π, 0) and (x0 , p0 ) = (0.6π, 0).

Note: some of the methods might not be stable (either |yn| - o or failure of the Newton method to converge) for all values Ni  and your code should be able to handle that (and you should remark on the cases where this occurred in your report).

Important:  like always, don’t forget to discuss how the results you observed relate to the theory. What do we know about the orders of the three methods and also have a look at your results from Q1.3 of this assignment.

Q 2.3.  [5/20]

The following extends Q2.2 but is a bit more involved, requiring the implementation of a new method without the help of the quizzes:

Introduction:

Consider an ODE of the form y\ = f(t, y) = My +F (t, y) where the assumption is that the main dynamics is given by the linear part My while the non-linearity F (t, y) = f(t, y) -My is only a small perturbation. We can obtain an ODE for z = exp(-Mt)y (using the matrix exponential):

z\ = -M exp(-Mt)y + exp(-Mt)y\ = exp(-Mt)F╱t, exp(Mt)z、=: G(t, z)

For this ODE we can compute an approximation zn  to the exact solution Y (tn) using for example a Runge-Kutta method. Then an approximation to the original solution Y (tn) is given by yn  := exp(Mtn)zn .

Note that if F = 0, i.e., f(t, y) = My then G(t, z) = 0 and zn  = z0  = y0  for all n and therefore yn = exp(Mtn)y0  which is the exact solution to y\ = My .

Scheme of that type are called exponentially fit  method or exponential  integrators.  The main property being that they solve linear problems exactly.   They are often used for Hamiltonian systems with a Hamiltonian close to linear.

Example:

Consider the Forward-Euler method for z\  = G(t, z), i.e., zn+1  = zn + hG(tn, zn).  Multi- plying with eMtn41 :

yn+1 = eMtn41 ╱zn + h exp(-Mtn)F╱tn , exp(Mtn)zn、∶

Note that eMtn41   = eMtn eMh  so

yn+1 = eMh ╱yn + hF╱tn, yn、∶ .

While the Backward-Euler method zn+1 = zn+ hG(tn+1, zn+1 for z\ = G(t, z) leads to the following method:

yn+1 = eMtn41zn+1 = eMh eMtnzn+heMtn41e-Mtn41F ╱tn+1 , exp(Mtn+1)zn+1、= eMh zn+hF╱tn+1, yn+1、

Important to note that, in both cases, we do not need to carry out the expensive compu- tation of the matrix exponential eMtn  in each time step but only need eMh which for fixed h we can compute in advance.

Method:

We will implement  an  exponentially fit  DIRK. These  are based on  a standard DIRK method with a Butcher tableau given by α, β, γ . The exponentially fit version is then

之                                                                          s

k之 = e-α_i丘 F ╱tu + α之 h, eα_i丘 ╱yu + h      β之d  kd、∶ ,   yu+1 = eα丘 ╱yu + h      γ之 k之 ∶ .

d=1                                                                          之=1

Again note that we can pre-compute the required matrix exponentials for a given step

size, i.e., (eα_i丘)之(s)=1 , ( ╱eα_i丘∶ -1 )之(s)=1 , and eα丘 .

Example:

Let’s look again at the Forward-Euler method with s = 1, α 1 = 0, γ1 = 1 and β11 = 0:

k1 = F ╱tu, yu、,

yu+1 = eα丘 ╱yu + hk1 ∶ = eα丘 ╱yu + hF╱tu, yu、∶ .

While the Backward-Euler method (s = 1, α 1 = 1, γ1 = 1, β11 = 1) leads to

k1 = e-α丘 F ╱tu + h, eα丘 ╱yu + hk1、∶ ,

yu+1 = eα丘 ╱yu + hk1 ∶