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.  [2/20]

Consider the ODE y\  = f (y) with initial condition y(0) = y0  and the following approxi-

mation

yn+1 = yn +  ╱ f ╱yn、+ f ╱yn + hf (yn)、、

We use the notation from the lecture: h is the step size and yn ≈ Y (tn) with tn = nh.

(1) Consider the method applied to the linear ODE y\  = λy with real valued λ < 0. For which values of λ is the method stable?

(2) Consider the method applied to the linear ODE y\  = λy with imaginary λ = iα . For which values of α is the method stable (i.e. the solution remains bounded)?

Q 1.2.  [6/20]

Consider the ODE y\ = f (y) with initial condition y(0) = y0  and the approximation from

the first question

yn+1 = yn +  ╱ f ╱yn、+ f ╱yn + hf (yn)、、

We use the same notation as in the lecture where h is the step size and yn  ≈ Y (tn) with tn = nh. We assume that f is Lipschitz and Y e C3 .

(1) Show that the truncation error satisfies τn = O(h3 ).

Hint: make sure you also Taylor expand the term f (y + hf (y)). Also make good use of the ODE!

(2) Using  τn   =  O(h3 )  and  Gromwall’s  lemma,  prove  that  the  method  converges quadratically:

en = O(h2 )       (for h small enough) .

Hint: use that f is Lipschitz to derive something like

|en+1| < R(hLf)|en| + τn

with a suitable polynomial R = R(µ).

Q 1.3.  [2/20]

Consider a scalar linear  ODE y\   =  λy  on  [0, T] with  λ  e  R fixed.   Assume that an approximation is given by

yn+1 = R(λh)yn

where R = R(µ) is a non-constant polynomial.

Prove that such an approximation is never unconditionally stable.

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 the following functions and test them individually.  Take a look at the quiz questions which also include some tests and more detail on the function signature to use. After the quiz closes you can use the provided solutions if you encountered problems implementing the required functions.

(1) Write a function that computes a single step of the forward Euler method

yn+1 = yn + hf (tn, yn)

given h, tn  e R and yn  e Rm  and a general vector valued function f which should be provided as a function handle to your function.

(The solution to this is available in the warmup quiz).

(2) Write a function that computes a single step of the method from Q 1.1 yn+1 = yn +  ╱ f (tn, yn) + f ╱tn + h, yn + hf (tn, yn)、、

given h, tn  e R and yn  e Rm  and a general vector valued function f which should be provided as a function handle to your function.

(3) Write a function evolve to solve a vector valued ODE of the form y\ (t) = f (t, y(t)) ,        t e (0, T) ,        y(t0 ) = y0  .

using an approximation of the form:

yn+1 = Φ(tn, yn; h) .

The code should be general with respect to the dimension of the ODE and the right hand side f .  The right hand side and method function Φ should both be parameters (function handles) to your function  evolve  together with t0 , y0 , T, N and use a fixed step size h =  .  The method should return the vectors (tn)n(N)=0 and (yn)n(N)=0 .

(4) Write a function computeEocs to compute the experimental order of convergence (EOC) for a given sequence (hi, Ei)i(r)=0 . So the function should return

log ╱ 、

eoci  = log ╱ 、  ,        i = 1, . . . , r

where Ei = Ehi   is the error of a method with step size hi .

Note: I explained the concept of EOC during the seminar in week 2 (recorder) and there is a short recorded lecture on this concept which will be uploaded during week 3. Basically eoci should be close to the convergence rate of the method tested (at least for h small enough).

Q 2.1.  [2/20]

Test your implementatioqn using the 2x2 ODE: y\ (t) = f (t, y(t)) and y(0) = y0  with

f (t, y1 , y2 ) = ╱ y2 (λy2- 2y1 ) \

for t e [0, T] with T = 10.

With λ = 1 and y0 = (2, -2)T  this has the exact solution

Y (t) =  ╱            、

．    !

Compute the maximum errors Ehi   and the EOCs for the sequence of time steps given by hi = N达(T)2i   for i = 0, . . . , 9 using N0 = 25.

Q 2.2.  [2/20]

Write a program to solve a vector valued ODE of the form

y\ (t) = f (t, y(t)) ,        t e (0, T) ,        y(0) = y0  ,

using the method from Q 1.1

yn+1 = yn +  ╱ f (tn, yn) + f ╱tn + h, yn + hf (tn, yn)、、  .

Use the same evolve  method as before but with a different parameter Φ .  Implementing this was part of the first quiz.

Make sure to optimized your implementation with respect to the number of calls to f . Repeat the experiments from Q 2.1.

Q 2.3.  [3/20]

Compare the results from Q 2.1 and Q 2.2 with respect to their efficiency. To check effi- ciency compare the errors of the two methods with respect to the number of f evaluations needed to compute the approximations y0 , . . . , yN . Make sure your implementation of the methods from Q 2.2 does not make any unnecessary evaluations of f !

Note:  you can but don’t need to change your code to keep track of the evaluations of f needed.  In this case it is sufficient to just explain the difference between the methods based on your implementation noting that for a given N , both methods call Φ N times. So, the only difference between the methods comes from the number of calls to f in your implementation of the different Φ functions.

Q 2.4.  [3/20]

Investigate the approximations given by both methods on the time interval [0 , 3] for the ODE y\ (t)  =  f (t, y(t)) using the  method  of manufactured  solutions  with the following exact solution

Y (t) =

sin(t)

et-

t <

t >