Description
1. Consider the ODE
dy(t)
dt = −3y(t) sin t, y(t = 0) = π
√
2
,
which has the exact solution y(t) = πe3(cos t−1)/
√
2 (you can verify that). Implement the methods
forward Euler and Heun’s for this ODE to test the error as a function of ∆t . In particular:
(a) Solve the ODE numerically using the forward Euler method:
y(tn+1) = y(tn) + ∆tf(tn, y(tn))
with t = [0 : ∆t : 5], where ∆t = 2−2
, 2
−3
, 2
−4
, . . . , 2
−8
. For each of these ∆t values calculate the
error E = mean(abs(ytrue −ynum)) of the numerical method. Plot log(∆t) on the x axis and log(E)
on the y axis. Using polyfit, find the slope of the best fit line through this data. This is the order
of the forward Euler method.
ANSWER: Save your last numerical solution (∆t = 2−8
) as a column vector in A1. Save the
error values in a row vector with seven components in A2. Save the slope of the line in A3.
(b) Solve the ODE numerically using Heun’s method:
y(tn+1) = y(tn) + ∆t
2
[f(tn, y(tn)) + f(tn + ∆t, y(tn) + ∆tf(tn, y(tn)))]
with t = [0 : ∆t : 5], where ∆t = 2−2
, 2
−3
, 2
−4
, . . . , 2
−8
. For each of these ∆t values, calculate
the error E = mean(abs(ytrue − ynum)) of the numerical method. Plot log(∆t) on the x axis and
log(E) on the y axis. Using polyfit, find the slope of the best fit line through this data. This is
the order of the Heun’s method.
ANSWER: Save your last numerical solution (∆t = 2−8
) as a column vector in A4. Save the
error values in a row vector with seven components in A5. Save the slope of the line in A6.
2. Consider the van der Pol oscillator
d
2y(t)
dt2
+ ǫ[y
2
(t) − 1]dy(t)
dt + y(t) = 0
with ǫ being a parameter.
(a) With ǫ = 0.1, solve the equation for t = [0 : 0.5 : 32] using ode45. The initial conditions are
y(t = 0) = √
3 and dy(t = 0)/dt = 1. Repeat this for ǫ = 1 and ǫ = 20.
ANSWER: Save the solutions y(t) for different ǫ as a matrix of 3 columns in A7.
(b) Using the time span t = [0, 32] (the step size for displaying the result is not specified), solve
the van der Pol’s equation with ode45. Use ǫ = 1 and the initial conditions y(t = 0) = 2 and
1
dy(t = 0)/dt = π
2
.
Below is an example on how to control the error tolerance TOL in ode45:
TOL = 1e-4;
options = odeset(‘AbsTol’,TOL,‘RelTol’,TOL);
[T,Y] = ode45(‘rhs’,tspan,y0,options);
Using the diff and mean commands on the vector T shown above, calculate the average stepsize t needed to solve the problem for each of the following tolerance values: 10−4
, 10−5
, . . . , 10−10
.
Plot log(∆t) on the x axis and log(TOL) on the y axis. Using polyfit, find the slope of the best fit
line through this data. This is the order of the local truncation error of ode45. Repeat this with
ode23 and ode113.
ANSWER: The slopes should be saved in variables A8 – A10 for ode45, ode23, and ode113 respectively.
3. To explore interaction between neurons, implement two Fitzhugh neurons coupled via linear
coupling:
dv1
dt = −v
3
1 + (1 + a1)v
2
1 − a1v1 − w1 + I + d12v2
dw1
dt = bv1 − cw1
dv2
dt = −v
3
2 + (1 + a2)v
2
2 − a2v2 − w2 + I + d21v1
dw2
dt = bv2 − cw2
with parameters a1 = 0.05, a2 = 0.25 , b = c = 0.01 and, I = 0.1. Start the simulations with the
initial condition of (v1(0), v2(0)) = (0.1, 0.1) and (w1(0), w2(0)) = (0, 0) and use the ode15s solver.
Set the interaction parameters such that d12 is negative and d21 is positive. What do you observe
from the different graphical representations of the solutions?
ANSWERS: Set the interaction parameters to 5 different values
(d12, d21): (0,0), (0,0.2), (-0.1,0.2), (-0.3,0.2), (-0.5,0.2). For each interaction value solve the system
for t = [0 : 0.5 : 100] and save the computed solution, (v1, v2, w1, w2), in a 201 × 4 matrix. Write
out the 5 different solutions, each of which corresponds to an interaction parameter, in A11 – A15.
