============ Part B – Main Part ========================

Fakult¨ at f¨ ur Mathematik und Informatik 29. Mai 2012 TU Bergakademie Freiberg

W. Queck/M. Helm

Numerical Analysis of Differential Equations Initial Value Problems (IV)

============ Part B – Main Part ========================

Exercise 1

Determine the solution of the following linear difference equations. For the identification of a par- ticular solution of the inhomogeneous difference equation try an “ansatz“ with the same structure as the right hand side.

a) x

=

(x

+ x

), n = 3, 4, . . ., Initial conditions: x

=

, x

=

b) x

− 2x

+ x

= 0

c) 2x

+ 3x

− 2x

= 0, n = 0, 1, . . ., Initial conditions: x

= 5, x

= 0

d) x

− 5x

+ 6x

= 2, n = 2, 3, . . ., Initial conditions: x

= 1, x

= −1 Exercise 2

Restricted 3-body problem

We consider a satellite in the gravity field of earth and moon. We assume that the movement 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 + µ

[(x + µ)

+ y

]

− µ x − µ [(x − µ)

+ y

]

, y

= y − 2x

− µ y

[(x + µ)

+ y

]

− µ y

[(x − µ)

+ y

]

.

Here, µ = 1/82.45 is the ratio 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:

y

(0) = 1.2, y

(0) = 0, y

(0) = 0, y

(0) = −1.04935750983031990726.

Therefore choose an appropriate end tend for the integration interval.

Initial Value Problems 2

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

– the tolerances RelTol and AbsTol for the stepsize control of ode45, – 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 horizon- tal 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 trajec- tory.

f) Solve the 3-body problem with a classical Runge-Kutta method of fourth order and constant stepsize h = (t

−t

)/n. How many integration steps n are needed to get a similar trajectory as with ode45?

Hint: 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.