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 two step LLM

yn+2 + ayn+1 + byn = hcf (yn+s)

with fixed constants a, b, c. Answer the following questions for s = 0, s = 1, and s = 2.

(i) Give conditions for a, b, c so that the method converges.

(ii) Find values for a, b, c so that the method converges and has consistency order of at least 2 (if such a method exists for a given value of s).

Q 1.2.  [3/20] Consider the non-linear oscillator given by

x\\ (t) + x(t) + ε ╱x2 (t) - 1←x\ (t) = 0 ,        x(0) = x0 , x\ (0) = 0

which is a non-dimensionalized form of the van der Pol equation. It can be shown that this problem has a solution that will (eventually) be periodic (the model has a so called limit cycle).  Since the problem is homogeneous we can start anywhere on the cycle, e.g., we can choose a point on the cycle with x\  = 0 (which should explain the initial conditions). The correct value for x0  to be on the cycle is unkown.

For the following you will need to be able to solve linear ODEs of the form y\\ + y = F (t) where the right hand side is of the form F (t) = F0 cos(ωt) (or F (t) = F0 sin(ωt)) for some ω . The general form of the solution is with constants of integration C1, C2 :

Y (t) :=    (

F (t) + C1 cos(t) + C2 sin(t)   ω  1 ,

t(t) + C1 cos(t) + C2 sin(t)        ω = 1 .

where \ = F . This should be known from MA133 or e.g. follow this link: Example 6.1.4 and Example 6.1.5.

Note that the problem is linear so a right hand side of the form A cos(ω1t) + B sin(ω2t) can be handled with the same formula, i.e., the solution is the sum of the solutions Y1, Y2 given above with right hand sides F (t) = A cos(ω1t) and F (t) = B sin(ω2t), respectively. The central thing to take into account is the difference in the solution for ω  1 and ω = 1 - the second case is know as secular solution.

(i) Using regular perturbation theory of the form x0 (t)+εx1 (t)+O(ε2 ) derive a system of ODEs for x0 , x1 .

(ii) Solve the ODE you obtain for x0 , x1  under the assumption that we are looking

only for periodic solutions (i.e. the limit cycle).

Q  1.3.  [4/20] Consider the problem from Q1.2.  Using an expansion up to O(ε3 ) and after substituting x0 , x1  into the ODE for x2  one obtains

x + x2 =  cos(t) + 2C sin(t) -  cos(3t) + 3C sin(3t) +  cos(5t)

with some free parameter C .

You can see that the right hand side contains two terms with cos(t) and sin(t) which lead to secular terms in x2 .  There is only one free constant C to choose which is not enough to eliminate both of these terms. So we can not obtain a periodic expansion above O(ε2 ). To construct a periodic expansion of higher order we need to allow the frequency of the oscillator to adapt to ε . To this end we define a stretched time variable τ = ωt where the frequency ω also has an expansion of the form ω = 1 + cω1 + O(ε2 ).

(i) Derive an ODE for the function y(τ ) = x().

(ii) Use the expansion for ω and y0 (t) + εy1 (t) + O(ε2 ) to derive a system of ODEs for y0 , y1  (which should depend on ω1 ).

(iii) Find periodic function y0 , y1  - what choice for ω 1  does this require?

Remark: Expanding both ω and y up to O(ε3 ) leads to an equation for y2  that looks so unpleasant that I didn’t want to force you to write it down.  After substituting y0 , y1 , ω 1 into the ODE for y2  one obtains

y + y2 = (4ω2 + ) cos(τ) + 2C sin(τ) -  cos(3τ) + 3C sin(3τ) +  cos(5τ)

with some constant C .  It is now possible to obtain a periodic solution y2  by choosing C, ω2  in such way that the secular terms are removed.

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.

Scipy initial value solver: Scipy provides a solver based on adaptive solvers as discussed in week 9. It quite easy to write an evolve type method using this solver:

from   scipy . integrate   import   solve _ ivp

def   evolveScipy ( method ,  f , Df ,  t0 , y0 ,  T ,  rtol ,   atol ):

res   =   solve _ ivp (f ,   [ t0 , T ] ,  y0 ,  method = method ,   rtol = rtol ,   atol = atol ) return   res . t ,  res . y . transpose ()

Note  that  the  Scipy  solver  returns  y  differently  from  our  setup which  can  be  solved by using transpose as shown above.   Have a look here to get an idea of the parame- ters to call this method https://docs.scipy.org/doc/scipy/reference/generated/ scipy.integrate.solve_ivp.html.  In this assignment we will use the method="RK45" which is similar to the embedded explicit Runge-Kutta method you will be implementing. The aim is to see if your code is competitive with the Scipy implementation.

Since these methods are adaptive they do not use a constant time step but one that it adapts to the solution. You can obtain the time steps used in the following way:

Think about ways to visualize the varying time steps by for example overlaying the loga- rithm of the time step with the actual solution on one graph or coloring the phase space

portrait according to the time step used - I would be interested to see what that looks like; so if you tried that and decided it is not so informative mention that in the report.

For these methods you do not chose a value of N (or h) to increase accuracy (or decrease computational cost) but you chose tolerances. Instead of having two different values for the tolerances (for the absolute and relative error), I suggest to combine the two by choosing one tolerance tol and then use for example rtol=tol,atol=tol*1e-2.

A lot of this assignment is about coming up with ways to visualise the outcome of sim- ulations and comparing different methods, which is as you hopefully figured out by now, a major challenge in computational mathematics. I’ve made some suggestions below you don’t need to follow all of them but some...

Q 2.1.  [7/20]

We will solve the van der Pol oscillator discussed in part 1 (and in the second part of the lecture in the context of Hamiltonians).

Start by having a look at the problem using the build-in method. You can do this question for a partial mark using only the build-in method.

Model:

Solve the van der Pol equation

x\\ - µ(1 - x )x2\ + x = 0

with µ = 200 which you of course will have to rewrite as a first order system.

We don’t know the exact solution but it is possible to obtain an estimate Θ for the period of the oscillation

Θ = ╱3 - 2 log(2) ←µ + 2πµ-    .

We will use T = 3.1Θ as final time so that we get three estimates for the period which we can average.

Here is a webpage discussing this problem including some Python code https://www. johndcook.com/blog/2019/12/22/van-der-pol. There is also a link discussing the cal- culation of the period - although you will have to adapt that for your own implementation.

Task:

Implement your own adaptive time stepping method - the idea is explained in the record- ing for week 9.  A short summary of the method is also given on the final page of this assignment. Compare your method with the build-in method.

Setup:

Compare your implementation with the build-in version by trying to obtain a relative error in the period of below 10-3 without too much computational cost. For each method, try to find a value for the tolerance tol that is as large as possible so that the methods produce solutions with a relative error in the period of less than 1e - 3.  Note that the build-in implementation is of course not identical to the one you implemented so you might need different tolerance to get a sufficiently accurate result. My suggestion is to automate finding the  “optimal” tolerance tol by starting with a fairly large value (not too large because the method might fail) and then decreasing this value by a fixed factor, e.g., 2, until the desired accuracy for the period is reached.

You will need some way to estimate the period Θ given a discrete approximation (tn, yn)n . You could for example find n1 , n2 , .., nk  where yn  crosses the x- or y-axis.  Then tni   - tni −2    is an approximation to the period for i  =  2, 3, 4, ..., k .   You could average these approximations to obtain an average discrete period Θh  produced by the approximation. First test your approach by seeing that you can get the build-in method to produce an approximation with |Θ - Θh|/|Θ| < 10-3 .

Discussion:

Once you obtained the “optimal“ tolerance for each method, discuss how the two method compare with respect to number of steps required.  Are the time step chosen over the time interval similar?  Also estimate the computational cost of both methods - you can use the Python time module for example or add a counter to the function f to find out the number of evaluations of the right hand side required (you can use a global variable which you increment each time f is called).

Final remark:

the build-in method has some way of finding a good value for the first h while we need to prescribe one - you could start with the same one that the build-in method uses. Also try to chose hmax large enough so that it does not influence the outcome - it should not be needed for this problem.

Q 2.2.  [3/20]

You can do this question only using the build-in method for partial marks if you didn’t manage to implement you own adaptive method.

Have a look at the problem described here:  https://www.johndcook.com/blog/2020/   02/08/arenstorf-orbit.   Try to reproduce the phase portrait included on this page   showing the stable periodic orbit which apparently was the basis for the Apollo mission.     This problem is also described in one of the books on the reading list (“Solving ordinary   differential equations“ by Hairer, Nørsett, and Wanner page 127). Use the more accurate   initial conditions given there y0 = (0.994, 0, 0, -2.00158510637908252240537862224). That   source also provides an approximation to the period Θ = 17 .0652165601579625588917206249. Chose T = 3Θ and try to find tolerances so that the final value is close to the initial con-   ditions, e.g., |uN - u0 |/|u0 | < 0.01. You might have to chose quite a small tolerance. Can   you manage to maintain a periodic orbit with T = 5Θ for example?