Assignment 3 Numerical Solution of ODEs

CS 370 Numerical Computation Assignment 3
Numerical Solution of ODEs
What you need to get
• https://www.overleaf.com/read/vztpkwrpmtcn: (optional) Overleaf template
• YOU_a3-q3.ipynb: Jupyter notebook for Q3
• YOU_a3-q4.ipynb: Jupyter notebook for Q4
• YOU_a3-q5.ipynb: Jupyter notebook for Q5
What you need to know
This assignment uses a very simple class called MyODESolution. Its definition is in the notebook for
question 3.
What to do
1. [2 marks] Prove that
x(t) = 3e
t − 2e
2t
y(t) = 3e
t − 4e
2t
is a solution to the IVP
x
0
(t) = y(t)
y
0
(t) = −2x(t) + 3y(t)
x(0) = 1 , y(0) = −1
Show your work.
2. [6 marks] Compute two Euler steps of size h = 0.1 for the initial value problem
dx
dt = xy − t + 2
dy
dt =
3x
y
+ 5 t x
x(1) = 2 , y(1) = −3
3. [10 marks] Complete the Python function MyODE,
sol = MyODE(de, tspan, y0, rtol=0.05, event=[])
See the function’s documentation for a description of the inputs. Your function should have the
following features:
(a) Modified Euler: Implement Modified Euler’s method.
c Jeff Orchard, Reinhold Burger 2019 v1.0 Page 1
CS 370 Numerical Computation Assignment 3
(b) Adaptive Time-Stepping: Do adaptive time-stepping so that your estimate of the local relative error is less than rtol for each step. If the error is less than or equal to rtol, increase
the time step by 20% for the next step. If the error is greater than rtol, then reduce your step
size by 50% and repeat the step.
(c) Events Function: Stop when an event function returns a negative value.
(d) Interpolate Terminal Event: If the integration ends because an event is triggered, your function should linearly interpolate between the last two points to find a more accurate estimate
for the time of the event. Then it should interpolate the state at that new, interpolated end
time. You may assume that the state and the event function both change linearly between the
last two time steps.
(Hint: Think similar triangles.)
4. [8 marks] A fox is chasing a rabbit on a field with an electric fence along one side. Suppose that
f(t) = (fx(t), fy(t)) is the location of the fox at time t (in metres), and (rx(t), ry(t)) is the location
of the rabbit (also a function of time). Suppose the fox always accelerates toward the rabbit. Let the
vector pointing from the fox to the rabbit be a(t, f) = (ax(t, f), ay(t, f)), as shown in the diagram.
Then, the initial value problem governing the fox’s
movements is,
f
00
x = A
ax(t,f)
ka(t,f)k
−

(f
0
x
, f0
y
)

 − S


f
0
x
(1)
f
00
y = A
ay(t,f)
ka(t,f)k
−

(f
0
x
, f0
y
)

 − S


f
0
y
(2)
fx(0) = 0 , fy(0) = 0 (3)
f
0
x
(0) = 0 , f0
y
(0) = 0 (4)
fox
y
x
electric fence
30 m
rabbit
a
where A and S are both positive scalar constants. Generally, A corresponds to the rate of acceleration of the fox, and S limits the fox’s maximum speed. The electric fence follows the line x = 30
(metres).
(a) Convert the system of differential equations to a system of first-order differential equations.
Show your work, and clearly state the new system. Place your solution in a markdown box
in the notebook using LATEX.
(b) Complete the function fox_dynamics in the notebook so that it computes the dynamics
function for your system of differential equations.
(c) The notebook includes the function rabbit(t), which gives the rabbit’s position as a function of time. Complete the event functions capture and fence to implement the following
two events:
i. capture: The first event is the fox catching the rabbit. You may assume that the fox
catches the rabbit if it comes within 10 cm of the rabbit.
ii. fence: The second event is the fox hitting the electric fence.
(d) Run some simulations of your system using parameter values in the range S ∈ [2, 2.5], and
A ∈ [4, 7]. You can either use your implementation of MyODE with a tolerance of 0.01, or you
can use solve_ivp. If you decide to use your version of MyODE, just copy the function (as
well as the definition of the class MyODESolution) into the notebook.
Find parameter values that result in each of the following outcomes:
c Jeff Orchard, Reinhold Burger 2019 v1.0 Page
CS 370 Numerical Computation Assignment 3
i. The fox hits the fence.
ii. The rabbit escapes.
iii. The fox catches the rabbit.
For each case, plot the rabbit’s trajectory and the fox’s trajectory on the same plot, along with
the fence.
5. [4 marks] Consider the IVP,
dy
dt =
4
1 + t
2
y(0) = 0
(a) Prove that y(t) = 4 arctan t is a solution to the IVP. (Don’t forget about the initial state.) Place
your solution in a markdown block in the notebook.
(b) Note that arctan 1 = π
4
. Thus, y(1) = π. Add code to your notebook to compute an estimate
of π by numerically solving the IVP. You may use your own MyODE function (by copying it),
or SciPy’s solve_ivp function.
(c) Calculate and display the relative error of your estimate to show that it is less than 10−8
.
What to submit
You must submit a series of PDF documents to Crowdmark. For each coding question, export your
jupyter notebook to PDF. When a proof or manual calculation is requested, you can typeset your solution
(eg. using LATEX or Word), or write your solution by hand. In either case, your solution should also be
submitted as a PDF. If you submit photos of handwritten work, it is your responsibility to ensure that
they are legible.
c Jeff Orchard, Reinhold Burger 2019 v1.0 Page 3