Due: 11:59:59pm Friday November 14th , 2014
AMATH 301 - Autumn 2014
Homework 4
VERY IMPORTANT: Remember to suppress all output and plotting before submitting your code.
Exercise 1: This exercise will have you studying how to numerically solve the following ODE:
y 0 (t) = cos(t)
y(0) = 0
for t starting at 0 and going to 5. First, solve this using Forward Euler with a timestep of ∆t = 0.1.
Store as a row vector where the j th column is your approximate solution at tj = (j − 1)∆t.
Answer: Should be written out as A1.dat.
Repeat this using a timestep of ∆t = 0.05 and then ∆t = 0.025. Plot your original solution, these
two solutions, and the exact solution of y(t) = sin(t), to verify that you have coded this correctly.
In fact, let’s do something better than ”eye-balling” it so that you don’t lose attempt(s) submitting
the wrong answer.
Compute
y − y|| =
q P the error between your approximations y˜(t) and the exact y(t) using the norm: ||˜
N
1
˜(tj ))2 . Denoting e0.1 , e0.05 , e0.025 as the errors for your approximations using
j=1 (y(tj ) − y
N
∆t = 0.1, ∆t = 0.05, ∆t = 0.025, respectively, compute the ratios e0.1 /e0.05 and e0.05 /e0.025 . Store
these ratios.
Answers: Should be written out as A2-A3.dat.
Note that since Forward Euler has a global error of O(∆t), the ratios should be roughly 0.1/0.05
and 0.05/0.025.
Now solve the ODE using RK2 with a timestep of ∆t = 0.1. Store your result as a row vector.
Answer: Should be written out as A4.dat.
Repeat this for ∆t = 0.05 and ∆t = 0.025. Obtain the corresponding errors using the same norm,
and then compute the ratios e0.1 /e0.05 and e0.05 /e0.025 . Again, store the ratios.
Answers: Should be written out as A5-A6.dat.
Note that since RK2 has a global error of O(∆t2 ), the ratios should be roughly (0.1/0.05)2 and
(0.05/0.025)2 .
1
AMATH 301 - Autumn 2014
Due: 11:59:59pm Friday November 14th , 2014
Homework 4
Exercise 2: Consider the damped pendulum equation (say the hinge needs oil, there’s friction in
the air, etc):
x00 (t) =
g
sin x(t) − δx0 (t),
L
where x(t) is the angle of displacement, g is the acceleration due to gravity, L is the length of the
pendulum, and δ is a damping coefficient. Note that this equation is non-linear. However, if we
assume that x 1, then we can approximate sin x(t) with x(t):
x00 (t) =
g
x(t) − δx0 (t)
L
Now that this equation is linear, it can be written in terms of a linear system by introducing
v(t) = x0 (t), i.e
d
dt
x
v
=
A11 A12
A21 A22
x
v
Find the matrix A if g = −10, L = 8, and δ = 3.
Answer: Should be written out as A7.dat.
Recall that Forward Euler can be written in terms of an iteration: xn+1 = (I + ∆tA)xn . Also recall
that we related the stability of Forward Euler with the spectral radius of M = I + ∆tA. Use these
concepts, and the matrix A above, to determine the value of ∆tmax so that Forward Euler is stable
for all ∆t < ∆tmax .
Answer: Should be written out as A8.dat.
Solve the linear system, with matrix A above, using Forward Euler from t = 0 until t = 50. For
initial conditions, use that [x(0), v(0)]T = [1, 0]T , and use ∆t = 0.95∆tmax for a timestep. Plot the
displacement x(t) to ensure that your solution is stable (albeit probably not very accurate because
the timestep is still quite large). Store the first 10 values of x(t), including the initial condition, in
a row vector.
Answer: Should be written out as A9.dat.
Now solve the same linear system, same initial conditions, same Forward Euler, only use ∆t =
1.05∆tmax for a timestep. Plot the displacement x(t) to verify the instability of Forward Euler at
a timestep just 5% above the limit. Store the first 10 values of x(t), including the initial condition,
in a row vector.
Answer: Should be written out as A10.dat.
You can verify that the scheme will be unstable even at a timestep just 1% above the limit.
2
AMATH 301 - Autumn 2014
Due: 11:59:59pm Friday November 14th , 2014
Homework 4
Exercise 3: For this exercise, you will be tracking a droplet as it moves around in a fluid vortex.
To do this, you’ll solve for the position z(t) = (x(t), y(t))T of particles along the droplet edge:
0
z (t) = f (t, z), where f (t, z) =
− sin2 (πx) sin(2πy) cos(πt/8)
sin2 (πy) sin(2πx) cos(πt/8)
The droplet itself is centered around [0.5, 0.75]T at t = 0 and has a radius of 0.15. Thus, to simulate
the droplet, we’ll track particles that start at zj (0) = [0.5 + 0.15 cos(θj ), 0.75 + 0.15 sin(θj )]T , where
j
θj = 99
2π, for j = 1, . . . , 100. First, solve for the position of these 100 marker particles at t = 2
using ode45. You may use the skeleton code markerparticle.m from the course website to get
started on this. If you do, make sure to upload all the required code to CSB.
Plot the position of the 100 marker particles at t = 2. Your particles should represent a deformed
droplet. Once you’ve verified this, store the positions of these particles at t = 2 as a 2x100 matrix
so that the first row is the x values, and the second row is the y values.
Answer: Should be written out as A11.dat.
Now, solve for the position of the 100 particles at t = 4. Plot the position of these particles. Your
droplet should have become more entrained in the fluid vortex. Once you’ve verified this, store the
positions at t = 4 as a 2x100 matrix.
Answer: Should be written out as A12.dat.
To continue, solve for the position of the 100 particles at t = 6. Plot the position of these particles.
The vortex has reversed direction at this point, and your droplet should be untraining. Once you’ve
verified this, store the positions at t = 6 as a 2x100 matrix.
Answer: Should be written out as A13.dat.
Finally, solve for the position of the 100 particles at t = 8. Plot the final position of these particles.
The droplet should look similar to how it started (albeit with quite a bit of numerical error at this
point). Once verified, store the positions at t = 8 as a 2x100 matrix.
Answer: Should be written out as A14.dat.
You just simulated this video https://www.youtube.com/watch?v=8V6kc0PQa14 with a computer.
Note that you can use odeset to set the relative error tolerance (RelTol ) lower in case you don’t
like the result at t = 8.
3