Description
1. (15 marks) Consider Newton’s equations of motion for a two-body problem specified by:
x
00(t) = −
x(t)
(x(t)
2 + y(t)
2)
3/2
; x(0) = 0.4; x
0
(0) = 0
y
00(t) = −
y(t)
(x(t)
2 + y(t)
2)
3/2
; y(0) = 0; y
0
(0) = 2.
As t ranges from 0 to 2π, (x(t), y(t)) defines an ellipse.
(a) Put this initial value problem into the standard form: z
0
(t) = f(t, z).
(b) Write the corresponding dynamics function in Matlab.
(c) On the course site you will find a Matlab function, y = ForwardEuler(f, t0,tf,N,z0),
which takes as input the system dynamics function f, the initial and final t values, t0
and tf , the number of steps, N, and the initial value, z0, and implements the forward
Euler method. Use this command to graph the resulting curves (x(t), y(t)) from t = 0
to t = 2π and verify that it indeed defines an ellipse.
2. ( 15 marks ) Consider the initial value problem
y
0
(t) = −2ty(t)
2
, y(0) = 1.
For this particular IVP it is easily checked that
y(t) = 1
1 + t
2
is the exact analytical solution. We will use this exact solution to determine the errors in
various numerical methods applied to solve the IVP.
(a) Write a Matlab function, y = MidpointRule(f,t0,tf,N,y0), which takes as input a
system dynamics function f, the initial and final t values, t0 and tf , the number of
steps, N, and the initial value, y0, and implements the midpoint Runge-Kutta method
(course notes on page 54). It should return the approximate solution as a vector y of
length N+1 such that
yi ≈ y(ti−1) for i = 1, 2, …, N + 1
where ti = t0 + ih, h =
tf −t0
N
, and y(ti−1) denotes the exact solution at t = ti−1.
(b) Write a Matlab script to apply the ForwardEuler function from the course web site to
solve the IVP to value tf = 4 with N = 40 time steps. Plot the analytical solution
(using solid lines) and the numerical solution (using the ‘.’ symbol for each point) on
the same graph, for 0 ≤ t ≤ 4.
1
(c) Repeat part (c) using the MidpointRule function.
(d) Calculate evidence illustrating the order of the forward Euler method as follows.
Let err(h) denote the error when computing with step size h. For a p-order method,
err(h) ≈ c · h
p
for some constant c, and so the effect of cutting the step size in half is
err(h/2)
err(h)
≈
1
2
p
i.e. the error is reduced by a factor of 2p
. Therefore, the order p of a given method can
be illustrated by calculating the above ratios.
For the forward Euler method, calculate a vector of successive computed values for the
solution at t = 4 based on cutting the step size in half a number of times. Specifically,
use N = 2i ×10 for i = 0, 1, …, 10. Calculate the errors and the successive ratios err(h/2)
err(h)
.
From your experimental evidence, what is the order of the forward Euler method?
(e) Repeat part (e) for the midpoint method. From your experimental evidence, what is the
order of the midpoint rule method?
3. (10 marks) Suppose we are solving the ordinary differential equation
y
0
(t) = f(t, y(t)) .
Determine the local truncation error of the method
yn+1 = 4yn − 3yn−1 − 2hf(tn−1, yn−1)
where yn is an approximation to y(tn), and tn = nh with h the stepsize.
4. ( 10 marks ) Consider the second-order Runge-Kutta method given by
yn+1 = yn +
h
4
f(tn, yn) + 3f
tn +
2h
3
, yn +
2
3
hf(tn, yn)
.
Analyse the stability of this method by applying the test equation
y
0
(t) = −λy(t), λ > 0.
What is your conclusion about stability limitations on h?
5. (25 marks) Consider a target moving in a two dimensional plane with location (xt
, yt). The
target is moving with constant speed St
, in direction given by the unit vector (cos θt
,sin θt),
where θt
is the angle of the velocity vector with respect to the x axis.
The target is pursued by a pursuer, with location (xp, yp). The pursuer is moving with
constant speed Sp, in direction given by the unit vector (cos θp,sin θp), where θp is the angle
of the velocity vector with respect to the x axis.
The pursuer attempts to hit the target, by turning as fast as possible towards the target. The
minimum turning radius of the pursuer is Rp, and the minimum turning radius of the target
is Rt
.
2
The trajectories of the pursuer and target can be described by the following system of ODEs,
where t is the time.
x
0
p
(t) = Sp cos θp(t), y0
p
(t) = Sp sin θp(t), θ0
p
(t) = Sp
θd(t) − θp
(|θd(t) − θp(t)| + D)
x
0
t
(t) = St cos θt(t), y0
t
(t) = St sin θt(t), θ0
t
(t) = St
Rt
ωt(t) .
where
(xp, yp) = x, y coordinates of pursuer
(xt
, yt) = x, y coordinates of target
θt = target moving in direction cos θtxˆ + sin θtyˆ
θp = pursuer moving in direction cos θpxˆ + sin θpyˆ
Sp = speed of pursuer
St = speed of target
Rt = minimum turning radius, target
θd = direction from pursuer to target (see Figure 1)
ωt = angular speed of target(fraction of maximum St)
D = turning damping factor, pursuer (1)
The damping factor D in equation (1) prevents the pursuer direction from changing too
rapidly when θd ‘ θp.
(xp
, yp
)
(xt, yt)
θd
Figure 1: Direction from pursuer to target.
3
Parameter Value
Initial Target Position (xt
, yt) = (5, 5)
Initial Pursuer Position (xp, yp) = (0, 0)
Initial Target Direction θt = 0.0
Initial Pursuer Direction θp = −π
Integration Interval [0, 20]
Rt
.25
Sp 2.0
St 1.0
Damping Factor D = 10−3
Hitting Distance dstop = 10−2
Angular Speed of Target ωt = 0.1
Error Tolerance tol = 10−6
Table 1: Data for pursuit problem.
Objective
The objective of this assignment is to develop Matlab code which solves the pursuit-target
problem, and determines the time (if it occurs) when the pursuer hits the target. The pursuer
is considered to have hit the target if the distance between the target and the pursuer is less
than dstop, i.e.
q
(xt − xp)
2 + (yt − yp)
2 ≤ dstop . (2)
The data for this problem is given in Table 1.
(a) (10 marks) Write Matlab code to solve the ODE system (1). Use an event function (see
course notes) to stop the simulation if the hitting criterion (2) is satisfied. In order to
determine the direction θd to the target (equation (1)), we have to be careful, since θp, θp
can be less than or larger than 2π (the pursuer/target could wind around a central point
several times). As well, we want to make sure the pursuer turns towards the target in
the shortest way possible. Use the supplied function direction.m (available on the course
web site) to compute θd.
Use the stiff solver ode15s. Use the following options to be passed to the ODE solver
(see the course notes).
options = odeset(’AbsTol’, tol,’RelTol’,tol,’MaxOrder’,5,’Stats’,’on’,…
’Events’,@event, ’Refine’,4);
Submit a hard copy of your code, a plot of the trajectories of the pursuer and target,
and the hitting time. The trajectory of the pursurer is given by the the parametric
curve (xp(t), yp(t)), and the trajectory of the target is given by the parametric curve
(xt(t), yt(t)) for t ∈ [0, T] where T is the stopping time or end of the integration interval
(whichever is first).
(b) (5 marks) Now, carry out some tests to see the effect of the error control tol on the
solution accuracy. Fill in Table 2. Explain what you see.
4
tol Hitting Time Number of Function
Evaluations
10−3
10−4
10−5
10−6
10−7
10−8
10−9
Table 2: Test of error control
(c) (5 marks) Now, try solving the problem again, except use ode45 as the solver. Use
tol = 10−4
, 10−5
, 10−6
, 10−7
. Compare with the solution using ode15s. Explain what
you see.
(d) (5 marks) Attempt to devise an evasion strategy for the target. The target is constrained so that St = 1.0, and the ωt ≤ 1. Some ideas: when the distance between the
target and the pursuer is less than a small multiple of Rt
, try turning the target in a
direction away from the pursuer. Increase the maximum integration interval to [0, 40].
Submit a pseudo-code description of your evasion strategy, a hard copy of your matlab
code, and a plot of the target and pursuer trajectories.
The five most interesting strategies (as determined by the TAs), will earn 3 bonus marks.
In the case of disputes, the TAs decisions are final.
5
