Lab 3 - Jahbriel Gann - MAT 275 Lab
Table of Contents
Exercise 1 .......................................................................................................................... 1
Exercise 2 .......................................................................................................................... 3
Exercise 3 .......................................................................................................................... 6
Exercise 4 .......................................................................................................................... 8
Exercise 5 .......................................................................................................................... 9
Introduction to Numerical Methods for Solving ODEs
Exercise 1
Part (a) Define ODE function f for ODE dy/dt = f(t,y) = y.
f=@(t,y)(y);
% Define vector t of time-values over the interval [0,2.5] to compute
% analytical solution vector.
t=linspace(0,2.5,100);
% Create vector of analytical solution values at corresponding t
values.
y=-2*exp(t);
[t100,y100]=euler(f,[0,2.5],-2, 100);
% Solve IVP numerically using forward Euler's method with 10
timesteps.
[t10,y10]=euler(f,[0,2.5],-2, 10);
% Solve IVP numerically using forward Euler's method with 100
timesteps.
[t100,y100]=euler(f,[0,2.5],-2, 100);
% Solve IVP numerically using forward Euler's method with 1000
timesteps.
[t1000,y1000]=euler(f,[0,2.5],-2, 1000);
% Solve IVP numerically using forward Euler's method with 10000
timesteps.
[t10000,y10000]=euler(f,[0,2.5],-2, 10000);
% Compute numerical solution error at t=2.5 for forward Euler with 10
% timesteps.
[t10,y10]=euler(f,[0,2.5],-2, 10);
y=2*exp(-2);
approx = y - y10;
% Compute numerical solution error at t=2.5 for forward Euler with 100
% timesteps.
[t100,y100]=euler(f,[0,2.5],-2, 100);
y=2*exp(-2);
appro = y - y100;
% Compute numerical solution error at t=2.5 for forward Euler with 700
% timesteps.
[t1000,y1000]=euler(f,[0,2.5],-2, 1000);
y=2*exp(-2);
appr = y - y1000;
1
Lab 3 - Jahbriel Gann - MAT 275 Lab
% Compute numerical solution error at t=2.5 for forward Euler with
7000
% timesteps.
[t10000,y10000]=euler(f,[0,2.5],-2, 10000);
y=2*exp(-2);
app = y - y10000;
% Compute ratio of errors between N=10 and N=100.
n100 = approx/appro (end);
% Compute ratio of errors between N=100 and N=1000.
n1000 = appro/appr (end);
% Compute ratio of errors between N=1000 and N=10000.
n10000 = appr/app (end);
clear
clc
f = @(t,y) y;
t = linspace(0,2.5,100); y = @(t) -2*exp(t);
yEnd = y(t(end));
ratio(1) = NaN;
N = [10,100,1000,10000];
for i = 1:numel(N)
[tN, yN] = euler(f,[0,2.5],-2,N(i));
yNend(i) = yN(end);
eN(i) = yEnd - yN(end);
if i>1
ratio(i) = eN(i-1)/eN(i);
end
end
%formating for tabling
format long g
N = N'; yNend = yNend'; eN = eN'; ratio = ratio';
TableOfResults = table(N,yNend,eN,ratio);
disp(TableOfResults)
N yNend eN ratio
_____ _________________ ____________________
________________
10 -18.6264514923096 -5.73853642909738
NaN
100 -23.6274327021244 -0.737555219282562
7.78048379168056
1000 -24.2890924486271 -0.0758954727798802
9.71803972315576
2
Lab 3 - Jahbriel Gann - MAT 275 Lab
10000 -24.3573763206299 -0.00761160077701817
9.97102646384616
Part (b). The consecutive errors decrease when the number of steps increases. As Euler's method step size
increases by a factor 10 the error decreases by a factor of 10 as well. This is shown by the ratio going
towards 10 meaning that the error is decreasing by about a factor of 10.
Part (c). Euler's method is an underestimate because the graph is concave up. We are making rectangles
from left endpoints. we are taking the smaller values which makes this a underestimate
Exercise 2
Part (a) Plot slopefield for dy/dt = -1.8y.
t = 0:.4:7; y = -9:1.5:6 ; % define grid of values in t and y
direction
[T,Y]=meshgrid(t,y); % creates 2d matrices of points in the ty-plane
dT = ones(size(T)); % dt=1 for all points
dY = -1.8*Y; % dy = -1.8*y; this is the ODE
quiver(T,Y,dT,dY) % draw arrows (t,y)->(t+dt, t+dy)
axis tight % adjust look
hold on
Part (b)
% Define vector t of time-values over the interval [0,7] to define
3
Lab 3 - Jahbriel Gann - MAT 275 Lab
% analytical solution vector.
t = linspace(0,7,100);
% Create vector of analytical solution values at corresponding t
values.
y = exp(-1.8*t);
% Plot analytical solution vector with slopefield from (a).
plot(t,y, 'k-', 'linewidth', 2);
Part (c)
% Define ODE function.
f=@(t,y)(5*y);
% Compute numerical solution to IVP with 7 timesteps using forward
Euler.
[t5,y5]=euler(f,[0,5],2, 5);
% Plot numerical solution with analytical solution from (b) and
slopefield
% from (a) using circles to distinguish between the approximated data
% (i.e., the numerical solution values) and actual (analytical)
solution.
plot(t5,y5, 'ro-', 'linewidth', 2)
hold off; % end plotting in this figure window
4
Lab 3 - Jahbriel Gann - MAT 275 Lab
The solution is not accurate because the dy = -1.8y dx so dy will change from positive to negative and
also keep growing to large values.
Part (d)
% Define new grid of t and y values at which to plot vectors for slope
% field.
figure
t = 0:.4:7; y = -.6:0.2:1.3;
% Plot slope field for dy/dt = -1.8y corresponding to the new grid.
[T,Y]=meshgrid(t,y); % creates 2d matrices of points in the ty-plane
dT = ones(size(T)); % dt=1 for all points
dY = -1.8*Y; % dy = -1.8*y; this is the ODE
quiver(T,Y,dT,dY) % draw arrows (t,y)->(t+dt, t+dy)
axis tight % adjust look
hold on
% Define vector t of time-values over the interval [0,5] to define
% analytical soution vector.
t = linspace(0,5,100);
% Create vector of analytical solution values at corresponding t
values.
y = 2*exp(-1.8*t);
% Plot analytical solution vector with slopefield from (a).
plot(t,y, 'k-', 'linewidth', 2);
5
Lab 3 - Jahbriel Gann - MAT 275 Lab
% Define ODE function.
f=@(t,y)(-1.8*y);
% Compute numerical solution to IVP with 10 timesteps using forward
Euler.
[t10,y10]=euler(f,[0,10],2, 20);
plot(t10,y10, 'ro-', 'linewidth', 2)
hold off;
Euler's method is closer to the actual line because the step size is smaller so it has large values from only
using 10 steps
Exercise 3
% Display contents of impeuler M-file.
type 'impeuler.m'
% Define ODE function for dy/dt = f(t,y) = y.
f=@(t,y)(1.5*y);
% Compute numerical solution to IVP with 10 timesteps using "improved
% Euler."
[t10,y10] = impeuler(f,[0,2.5],-2,10);
[t10,y10]
function [t,y] = impeuler(f,tspan,y0,N)
6
Lab 3 - Jahbriel Gann - MAT 275 Lab
% Solves the IVP y' = f(t,y), y(t0) = y0 in the time interval tspan =
[t0,tf]
% using Improved Euler's method with N time steps.
% Input:
% f = name of inline function or function M-file that evaluates
the ODE
% (if not an inline function, use: impeuler(@f,tspan,y0,N))
% For a system, the f must be given as column vector.
% tspan = [t0, tf] where t0 = initial time value and tf = final time
value
% y0 = initial value of the dependent variable. If solving a
system,
% initial conditions must be given as a vector.
% N = number of steps used.
% Output:
% t = vector of time values where the solution was computed
% y = vector of computed solution values.
m = length(y0);
t0 = tspan(1);
tf = tspan(2);
h = (tf-t0)/N; % evaluate the time step size
t = linspace(t0,tf,N+1); % create the vector of t values
y = zeros(m,N+1); % allocate memory for the output y
y(:,1) = y0'; % set initial condition
for n=1:N
y(:,n+1) = y(:,n) + (h/2)*(f(t(n),y(:,n))+(f(t(n+1), y(:,n) +
h*f(t(n),y(:,n)))));% Implement Improved Euler's method
end
t = t'; y = y'; % change t and y from row to column vectors
end
ans =
0 -2
0.25 -2.890625
0.5 -4.1778564453125
0.75 -6.03830814361572
1 -8.7272422388196
1.25 -12.613592298294
1.5 -18.230582618628
1.75 -26.3488889409857
2 -38.0823785475185
2.25 -55.0409377444603
2.5 -79.5513553337902
How does your output compare to that in the protocol? The improved euler 100 step was less accurate than
the 10 step euler's method because triangles are more accurate than rectangles
7
Lab 3 - Jahbriel Gann - MAT 275 Lab
Exercise 4
Part (a) Define ODE function f for ODE dy/dt = f(t,y) = y.
f=@(t,y)(y);
% Define vector t of time-values over the interval [0,2.5] to compute
% analytical solution vector.
t=linspace(0,2.5,100);
% Create vector of analytical solution values at corresponding t
values.
y=-2*exp(t);
[t100,y100]=impeuler(f,[0,2.5],-2, 100); %Use Improved Euler's method
% Solve IVP numerically using forward Improved Euler's method with 10
timesteps.
[t10,y10]=impeuler(f,[0,2.5],-2, 10);
% Solve IVP numerically using forward Improved Euler's method with 100
timesteps.
[t100,y100]=impeuler(f,[0,2.5],-2, 100);
% Solve IVP numerically using forward Improved Euler's method with
1000 timesteps.
[t1000,y1000]=impeuler(f,[0,2.5],-2, 1000);
% Solve IVP numerically using forward Improved Euler's method with
10000 timesteps.
[t10000,y10000]=impeuler(f,[0,2.5],-2, 10000);
% Compute numerical solution error at t=2.5 for forward Improved Euler
% with 10 timesteps.
y=2*exp(-2);
approx = y - y10 (end);
% Compute numerical solution error at t=2.5 for forward Improved Euler
with 100
% timesteps.
y=2*exp(-2);
appro = y - y100 (end);
% Compute numerical solution error at t=2.5 for forward Improved Euler
with 1000
% timesteps.
y=2*exp(-2);
appr = y - y1000 (end);
% Compute numerical solution error at t=2.5 for forward Improved Euler
with 10000
% timesteps.
y=2*exp(-2);
app = y - y10000 (end);
% Compute ratio of errors between N=10 and N=100.
n100 = approx/appro (end);
% Compute ratio of errors between N=100 and N=1000.
n1000 = appro/appr (end);
% Compute ratio of errors between N=1000 and N=10000.
8
Lab 3 - Jahbriel Gann - MAT 275 Lab
n10000 = appr/app (end);
clear
clc
f = @(t,y) y;
t = linspace(0,2.5,100); y = @(t) -2*exp(t);
yEnd = y(t(end));
ratio(1) = NaN;
N = [10,100,1000,10000];
for i = 1:numel(N)
[tN, yN] = impeuler(f,[0,2.5],-2,N(i));
yNend(i) = yN(end);
eN(i) = yEnd - yN(end);
if i>1
ratio(i) = eN(i-1)/eN(i);
end
end
%formating for tabling
format long g
N = N'; yNend = yNend'; eN = eN'; ratio = ratio';
TableOfResults = table(N,yNend,eN,ratio);
disp(TableOfResults)
N yNend eN ratio
_____ _________________ _____________________
________________
10 -23.8434326685287 -0.521555252878283
NaN
100 -24.358761448423 -0.00622647298394696
83.7641557624922
1000 -24.3649245898505 -6.33315563973724e-05
98.3154897517297
10000 -24.3649872870211 -6.3438584518849e-07
99.8312886041069
Part (b). As the improved euler's step size increases by a factor of 10 the error decreases by a factor of 100,
this can be shown by the the ratio getting closer and closer to 100 which means that the error is decreasing
by a factor of about 100
Exercise 5
Part (a)
9
Lab 3 - Jahbriel Gann - MAT 275 Lab
figure
% Plot slopefield for dy/dt = -1.8y.
t = 0:.4:7; y = -9:1.5:6 ; % define grid of values in t and y
direction
[T,Y]=meshgrid(t,y); % creates 2d matrices of points in the ty-plane
dT = ones(size(T)); % dt=1 for all points
dY = -1.8*Y; % dy = -1.8*y; this is the ODE
quiver(T,Y,dT,dY) % draw arrows (t,y)->(t+dt, t+dy)
axis tight % adjust look
hold on
Part (b)
% Define vector t of time-values over the interval [0,3] to define
% analytical solution vector.
t = linspace(0,3,100);
% Create vector of analytical solution values at corresponding t
values.
y = exp(-1.8*t);
plot(t,y, 'k-', 'linewidth', 2);
10
Lab 3 - Jahbriel Gann - MAT 275 Lab
Part (c)
% Define ODE function.
f=@(t,y)(y);
% Compute numerical solution to IVP with 5 timesteps using forward
Improved Euler.
[t5,y5]=impeuler(f,[0,3],2, 5);
plot(t5,y5, 'ro-', 'linewidth', 2)
hold off; % end plotting in this figure window
11
Lab 3 - Jahbriel Gann - MAT 275 Lab
The neumerical solution is not accurate because the dy = -3y dx. The dy will keep changing positive to
negative and keep growing to larger values due to the step size not being large enough to keep the values
from growing exponentially.
Part (d)
% Define new grid of t and y values at which to plot vectors for slope
% field.
figure %New figure
t = 0:.4:7; y = -.6:.2:1.3;
% Plot slope field for dy/dt = -1.8y corresponding to the new grid.
[T,Y]=meshgrid(t,y); % creates 2d matrices of points in the ty-plane
dT = ones(size(T)); % dt=1 for all points
dY = -1.8*Y; % dy = -1.8*y; this is the ODE
quiver(T,Y,dT,dY) % draw arrows (t,y)->(t+dt, t+dy)
axis tight % adjust look
hold on
% Define vector t of time-values over the interval [0,3] to define
% analytical soution vector.
t = linspace(0,3,100);
% Create vector of analytical solution values at corresponding t
values.
y = exp(-1.8*t);
% Plot analytical solution vector with slopefield from (a).
plot(t,y, 'k-', 'linewidth', 2);
% Define ODE function.
12
Lab 3 - Jahbriel Gann - MAT 275 Lab
f=@(t,y)(-1.8*y);
% Compute numerical solution to IVP with 20 timesteps using forward
Improved Euler.
[t100,y100]=impeuler(f,[0,3],2, 100);
plot(t100,y100, 'ro-', 'linewidth', 2)
hold off;
Improved Euler's method more accurate to the actual line this time due to the step size being much smaller
so it brings in the large values from using 10 step
13
