Consider the so-called R¨ossler system, given in the form of a system of three ordinary differential equations as

  x˙   =   −(y + z),

〈   y˙   =   x +  y,                                                                         (1)

(  z˙   =    5(1)  + z (x − µ).

where (x,y,z) ∈ R3  and µ ∈ R is a parameter. For this assignment, we consider µ = 4 fixed, for which the R¨ossler system has a (not necessarily unique) periodic orbit, denoted Γ .

Q1.  Adaptive stepsize integration of an IVP

The point u(0)  := (0, 4.5, 2) ∈ R3  lies rather close to this periodic orbit Γ, and someone told you that the period TΓ  ≈ 6.  With this question, you will explore the accuracy of numerical simulation with an adaptive stepsize integrator from the MATLAB suite.

(a) Use ode45 to approximate the solution of system (1) that starts at u(0)  over the time

interval [0, 6]. Set the options for odeset as suggested in Exercise sheet 10:

opts  =  odeset( 'reltol ' ,  1e−6,  'abstol ' ,  1e−6,  'refine ' ,  1);

. How many integration steps does ode45 take?

. Plot the x-coordinate of the solution version time.

. In a separate figure, plot the integration step sizes as follows:

plot(cumsum(diff(t)),  diff(t),  ' . − ' )

Here, t is the time array provided by the output from ode45; you may wish to change the default values for markersize and linewidth to your liking.

(b)  Suppose that, instead of the order-4 method ode45, you use an order-2 method or an

order-8 method.   How many integration steps do you expect these oder-2 and order-8 methods to take? Make sure to include arguments for your answer.

(c) Repeat part (a) using ode23, and again using ode89, instead of ode45; do keep the same options settings and include the plots in the same figures used for part (a).

. Are the number of integration steps used by ode23 and ode89 as expected?  Explain

your answer, where you should also provide a reason for why they are different or not from expected.

. For each of ode23, ode45, and ode89 provide a brief discussion (150 words max) ex-

plaining how and why the steps-size adapts during the first 10 integration steps.

. Observe a sudden drop in the integration step size near t = 6 for all three methods.

What happens here?

[10 marks]

Q2.  BVP with collocation

The periodic orbit Γ of system (1) with µ = 4 has also been approximated via pseudo-arclength continuation, where time was rescaled by the period TΓ  of Γ . The file roessler mu4 .dat on Canvas contains the approximation for Γ, which you can load into MATLAB with the command

>>  M  =  load( 'roessler mu4 .dat ' );

This command saves the data as a matrix M with the first column listing the (scaled) times t, and the following three columns the x-, y-, and z-coordinates.  The goal of this question is to find different approximations of the period TΓ  and of the point γ0  ∈ Γ at which Γ crosses x = 0 in the direction of decreasing x. Where necessary, round your answers to 6 decimal places.

The data in the file roessler mu4 .dat was obtained using degree-m orthogonal collocation over N mesh intervals. The degree of the piecewise-polynomial approximation is m = 4.

(a) How many mesh intervals are there?

. Plot the periodic orbit in (x,y,z)-space choosing a suitable view point; . Highlight the end points of the mesh intervals;

. Highlight the first mesh interval;

. Identify j and highlight in the same graph the mesh interval [tj −1,tj ] that contains the

point γ0 .

(b) Determine the values of the time interval [tj −1,tj ] representing the jth  mesh interval that contains γ0 , and find the corresponding (vector-valued) polynomial P(北,y,z)(t) = a0 + a1 t+ a2 t2 + a3 t3 + a4 t4  (round your coefficients to 4d.p.).

(i) What are tj −1  and tj , so what is the mesh size tj −   − tj −  , i = 1, . . . ,m used for the (extended) interior mesh points?

(ii) Plot the x-coordinates of the data points for this mesh interval versus t; Make sure

the data points are clearly highlighted and connect them with straight lines.

(iii) In the same figure, plot the graph of the x-coordinate of P(北,y,z)(t) for tj −1  ≤ t ≤ tj ;

here, you should use a relatively large number of t-values to represent the interval [tj −1,tj ], so the graph looks smooth.

(c)  Consider again the mesh interval [tj −1,tj ] from part (b).

(i) Determine the location of the Gauss collocation points zj,i  in [tj −1,tj ] and evaluate the associated (x,y,z)-coordinates using P(北,y,z)(t) from part (b).

(ii) For each of these zj,i verify that the derivatives of P(北,y,z)(t) match those of the vector

field (to a certain accuracy).

. Note that your vector field is scaled, so you can compare the direction of these

vectors, but not their lengths. Explain how you verified the match.

. The scaling factor applied to the vector field is precisely the period TΓ  of the Γ .

Determine TΓ  and explain to which accuracy you trust your answer.

[15 marks]