EXERCISES 1
%Part (a) Define ODE function f for ODE dy/dt = f(t,y) = 2y.
f = inline('2*y','t','y');
% Define vector t of time-values over the interval [0,0.5] to compute
% analytical solution vector.
t = linspace(0,0.5,100);
% Create vector of analytical solution values at corresponding t
values.
y = 3*exp(2*t);
% Solve IVP numerically using forward Euler's method with 5 timesteps.
[t5,y5] = euler(f,[0,0.5],3,5);
% Solve IVP numerically using forward Euler's method with 50
timesteps.
[t50,y50] = euler(f,[0,0.5],3,50);
% Solve IVP numerically using forward Euler's method with 500
timesteps.
[t500,y500] = euler(f,[0,0.5],3,500);
% Solve IVP numerically using forward Euler's method with 5000
timesteps.
[t5000,y5000] = euler(f,[0,0.5],3,5000);
% Compute numerical solution error at t=0.5 for forward Euler with 5
% timesteps.
e5 = y(end) - y5(end)
% Here, we define error as |analytical solution value - numerical
solution
% value|.
% Compute numerical solution error at t=0.5 for forward Euler with 50
% timesteps.
e50 = y(end) - y50(end)
% Compute numerical solution error at t=0.5 for forward Euler with 500
% timesteps.
e500 = y(end) - y500(end)
% Display table of errors (and ratios of consecutive errors) for the
% numerical solution at t=0.5 (the last element in the solution
vector).
disp('----------------------------------------------')
disp('| N | approximation | error | ratio |')
disp('|------|-----------------|----------|--------|')
disp('| 5 | 7.4650 | 0.6899 | N/A |')
1
disp('| 50 | 8.0748 | 0.0801 | 8.6148 |')
disp('| 500 | 8.1467 | 0.0081 | 9.8381 |')
disp('| 5000 | 8.1540 |8.1534e-04| 9.9835 |')
disp('----------------------------------------------')
e5 =
0.6899
e50 =
0.0801
e500 =
0.0081
----------------------------------------------
| N | approximation | error | ratio |
|------|-----------------|----------|--------|
| 5 | 7.4650 | 0.6899 | N/A |
| 50 | 8.0748 | 0.0801 | 8.6148 |
| 500 | 8.1467 | 0.0081 | 9.8381 |
| 5000 | 8.1540 |8.1534e-04| 9.9835 |
----------------------------------------------
EXERCISES 2
t = 0:.45:10; y = -30:6:42 ; % 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 = -2*Y; % dy = -2*y; this is the ODE
quiver(T,Y,dT,dY) % draw arrows (t,y)->(t+dt, t+dy)
axis tight % adjust look
hold on
f=inline('-2*y','t','y');
t=linspace(0,10,200); y=3*exp(-2*t); % defines the exact solution of
the ODE
[t200,y200]=euler(f,[0,10],3, 200); % solves the ODE using Euler with
200 steps
[t8,y8]=euler(f,[0,10],3, 8); % solves the ODE using Euler with 8
steps
plot(t200,y200,'ko-','linewidth',2)
plot(t8,y8,'ro-','linewidth',2)
figure
t = 0:.4:10; y = -1:0.4:3; % 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
2
dY = -2*Y; % dy = -2*y; this is the ODE
quiver(T,Y,dT,dY) % draw arrows (t,y)->(t+dt, t+dy)
axis tight % adjust look
hold on
[t16,y16]=euler(f,[0,10],3, 16); % solves the ODE using Euler with 16
steps
plot(t200,y200,'ko-','linewidth',2)
plot(t16,y16,'ro-','linewidth',2)
%Part (b). With the stepsize decreasing by the factor of k (k=10), and
the rotios are about 10, the error will
%change and decrease by the same factor of k (k=10).
%Part (c). the solution is concave up, and it hsould be understimate.
3
EXERCISES 3
f=inline('2*y','t','y');
t=linspace(0,.5,100); y=3*exp(2*t); % defines the exact solution of
the ODE
[t5,y5] = impeuler(f,[0,.5],3,5); % use @f if defined in separate
function
[t5,y5];
[~,~] = impeuler(f,[0,.5],3,50); %solves the ODE using Euler with 50
steps
%compare with the improved Euler approximation with 5 steps is already
more accurate than the Euler approximation with 50 steps!
EXERCISES 4
f=inline('2*y','t','y');
t=linspace(0,.5,100); y=3*exp(2*t); % defines the exact solution of
the ODE
[t5,y5] = impeuler(f,[0,.5],3,5); % use @f if defined in separate
function
[t50,y50] = impeuler(f,[0,.5],3,50); % use @f if defined in separate
function
[t500,y500] = impeuler(f,[0,.5],3,500); % use @f if defined in
separate function
4
[t5000,y5000] = impeuler(f,[0,.5],3,5000); % use @f if defined in
separate function
e5 = y(end) - y5(end); % error with N = 5
e50 = y(end) - y50(end); % error with N = 50
e500 = y(end) - y500(end); % error with N = 500
e5000 = y(end) - y5000(end); % error with N = 5000
ratio1 = e5/e50 ;
ratio2 = e50/e500;
ratio3=e500/e5000;
% Display table of errors (and ratios of consecutive errors) for the
% numerical solution at t=0.5 (the last element in the solution
vector).
disp('----------------------------------------------')
disp('| N | approximation | error | ratio |')
disp('|------|-----------------|----------|--------|')
disp('| 5 | 8.1081 | 0.0467 | N/A |')
disp('| 50 | 8.1543 |5.3555e-04 | 87.2394|')
disp('| 500 | 8.1548 |5.4284e-06| 98.6567|')
disp('| 5000 | 8.1548 |5.4357e-08| 99.8650|')
disp('----------------------------------------------')
----------------------------------------------
| N | approximation | error | ratio |
|------|-----------------|----------|--------|
| 5 | 8.1081 | 0.0467 | N/A |
| 50 | 8.1543 |5.3555e-04 | 87.2394|
| 500 | 8.1548 |5.4284e-06| 98.6567|
| 5000 | 8.1548 |5.4357e-08| 99.8650|
----------------------------------------------
EXERCISES 5
Part (a) Plot slopefield for dy/dt = -2y.
t = 0:.45:10; y = -30:6:42 ; % 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 = -2*Y; % dy = -2*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,10] to define
% analytical soution vector.
t = linspace(0,10,200);
% Create vector of analytical solution values at corresponding t
values.
y = 3*exp(-2*t);
% Plot analytical solution vector with slopefield from (a).
plot(t,y,'k-','linewidth',2);
% Part (c)
% Define ODE function.
5
f = inline('-2*y','t','y');
% Compute numerical solution to IVP with 8 timesteps using impEuler.
[t8,y8]= impeuler(f,[0,10],3,8);
% 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(t8,y8,'ro-','linewidth',2);
ylim([-30,42]);
hold off; % end plotting in this figure window
% Part (d)
% Define new grid of t and y values at which to plot vectors for slope
% field.
t = 0:.4:10; y = -1:0.4:3;
% Plot slope field for dy/dt = -2y 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 = -2*Y; % dy = -2*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,10] to define
% analytical soution vector.
t = linspace(0,10,200);
% Create vector of analytical solution values at corresponding t
values.
y = 3*exp(-2*t);
% Define ODE function.
f = inline('-2*y','t','y');
% Compute numerical solution to IVP with 16 timesteps using impEuler.
[t16,y16]= impeuler(f,[0,10],3,16);
plot(t,y,'k-','linewidth',2);
plot(t16,y16,'ro-','linewidth',2);
ylim([-1,3]);
hold off; % end plotting in this figure window
6
