Lab 3 - Nathan Florean - MAT 275 Lab
Table of Contents
Exercise 1 .......................................................................................................................... 1
Exercise 2 .......................................................................................................................... 2
Exercise 3 .......................................................................................................................... 5
Exercise 4 .......................................................................................................................... 6
Exercise 5 .......................................................................................................................... 7
Introduction to Numerical Methods for Solving ODEs
Exercise 1
Part (a)
f=@(t,y)(3*y);
clear t y t5 t50 y50 t500 y500 t5000 y5000
t=linspace(0,.5,100); y=2*exp(3*t); %evalute exact solution
[t5,y5] = euler (f,[0,.5],2,5); %evalute solution N=5
[t50,y50] = euler(f,[0,.5],2, 50); %evalute solution N=50
[t500, y500] = euler(f,[0,.5],2,500); %evalute solution N=500
[t5000, y5000] = euler(f,[0,.5],2,5000); %evalute solution N=5000
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 = 'N/A';
ratio2 = e5/e50;
ratio3 = e50/e500;
ratio4 = e500/e5000;
disp('----------------------------------------------')
disp('| N | approximation | error | ratio |')
disp('|------|-----------------|----------|--------|')
disp('| 5 | 7.4259 | 1.5375 | N/A |')
disp('| 50 | 8.7678 | 0.1956 | 7.8619 |')
disp('| 500 | 8.9433 | 0.0201 | 9.7273 |')
disp('| 5000 | 8.9614 | 0.0020 | 9.9720 |')
disp('----------------------------------------------')
----------------------------------------------
| N | approximation | error | ratio |
|------|-----------------|----------|--------|
| 5 | 7.4259 | 1.5375 | N/A |
| 50 | 8.7678 | 0.1956 | 7.8619 |
1
Lab 3 - Nathan Flo-
rean - MAT 275 Lab
| 500 | 8.9433 | 0.0201 | 9.7273 |
| 5000 | 8.9614 | 0.0020 | 9.9720 |
----------------------------------------------
Part (b). The error decreases as the stepsize decreases, and it does so portortionally. As the stepsize de-
creases by a factor of 10, the error decreases proporttionally by what it seems to be approaching by 10,
which is the same factor by which the stepsize decreases.
This confirms Eulor's method is of order h. Everytime the stepsize is decrease by a specifiec factor, the
error is also reduced by approximately the same factor.
Part (c). The exact graph of the Initial Value Problem y0 = 3y = f(t,y) that follows the equation y= 2e^-3t,
which is an exponential graph that is positive. This mean that a tangent line placed at any point on the
graph is always going extend underneath the graph in both directions. Since Euler's Method calculates the
values of the tangent line at time tn, the approximation will always underestimate in this case.
Exercise 2
Part (a) Plot slopefield for dy/dt = -3y.
t = 0:.3:5; y = -20:2:25 ;
[T, Y]= meshgrid (t,y);
dT = ones(size(T));
dY = -3*Y;
quiver(T,Y,dT,dY)
axis tight
hold on
2
Lab 3 - Nathan Flo-
rean - MAT 275 Lab
Part (b)
t= linspace(0,5,100);
y = 2*exp(-3*t);
plot(t,y,'k-','linewidth',2);
Part (c)
f=@(t,y)(-3*y);
[t6,y6] = euler(f,[0,5],2,6)
plot(t6,y6,'ro-','linewidth',2);
hold off; % end plotting in this figure window
t6 =
0
0.8333
1.6667
2.5000
3.3333
4.1667
5.0000
y6 =
3
Lab 3 - Nathan Flo-
rean - MAT 275 Lab
2.0000
-3.0000
4.5000
-6.7500
10.1250
-15.1875
22.7813
The reason that the numerical solution is inaccurate is because there are not enough points and the step
size is to large.
Part (d)
figure('Name','euler N=12');
% Plot slopefield for dy/dt =-3y
t = 0:0.3:5; y = -1:0.4:2;
[T,Y]=meshgrid(t,y);
dT = ones(size(T));
dY = -3*Y;
quiver(T,Y,dT,dY);
axis tight
hold on
t = linspace(0,5,100);
4
Lab 3 - Nathan Flo-
rean - MAT 275 Lab
y = 2*exp(-3*t);
plot(t,y,'k-','linewidth',2);
t = 0:0.3:5;
y = -1:0.4:2;
f=@(t,y)(-3*y);
[t12,y12] = euler(f,[0,5],2,12);
plot(t12,y12,'ro-','linewidth',2);
hold off
The graph is more accurate with the increase in the step size.
Exercise 3
% Display contents of impeuler M-file.
type 'impeuler.m'
f=@(t,y)(3*y);
clear t5 y5;
[t5,y5] = impeuler(f,[0,.5],2,5); %use @f if defined in sperate
function
[t5,y5]
5
Lab 3 - Nathan Flo-
rean - MAT 275 Lab
function [t,y] = impeuler(f,tspan,y0,N)
% Solves the IVP y' = f(t,y), y(t0) = y0 in the time interval tspan =
[t0,tf]
% using 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
% implement Improved Euler's method
y(:,n+1) = y(:,n) + (h/2)*(f(t(n),y(:,n)) + f(t(n)+h,y(:,n) +
h*f(t(n),y(:,n))));
end
t = t'; y = y'; % change t and y from row to column vectors
end
ans =
0 2.0000
0.1000 2.6900
0.2000 3.6181
0.3000 4.8663
0.4000 6.5451
0.5000 8.8032
The ouput matches the protocol perfectly
Exercise 4
% Part (a)
f = @(t,y) (3*y);
clear t y; clear t5 y5; clear t50 y50; clear t500 y500;
clear t5000 y5000;
6
Lab 3 - Nathan Flo-
rean - MAT 275 Lab
t=linspace(0,.5,100); y=2*exp(3*t); %evalute exact solution
[t5,y5] = impeuler(f,[0,.5],2,5); %evalute solution N=5
[t50,y50] = impeuler(f,[0,.5],2, 50); %evalute solution N=50
[t500, y500] = impeuler(f,[0,.5],2,500); %evalute solution N=500
[t5000, y5000] = impeuler(f,[0,.5],2,5000); %evalute solution N=5000
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 = 'N/A';
ratio2 = e5/e50;
ratio3 = e50/e500;
ratio4 = e500/e5000;
disp('----------------------------------------------')
disp('| N | approximation | error | ratio |')
disp('|------|-----------------|----------|--------|')
disp('| 5 | 8.8032 | 0.1602 | N/A |')
disp('| 50 | 8.9614 | 0.0020 | 81.2293 |')
disp('| 500 | 8.9634 | 2.0122 E-5| 97.9866 |')
disp('| 5000 | 8.9634 | 2.0163 E-7| 99.7976 |')
disp('----------------------------------------------')
----------------------------------------------
| N | approximation | error | ratio |
|------|-----------------|----------|--------|
| 5 | 8.8032 | 0.1602 | N/A |
| 50 | 8.9614 | 0.0020 | 81.2293 |
| 500 | 8.9634 | 2.0122 E-5| 97.9866 |
| 5000 | 8.9634 | 2.0163 E-7| 99.7976 |
----------------------------------------------
Part (b) The ratio of error is approximately 100, which is 10 times the factor of N. This provers that the
Improved Euler's method s of order h^2.
Exercise 5
% Part (a)
% Plot slopefield for dy/dt = -3y.
t = 0:.3:5; y = -20:2:25 ;
[T,Y]=meshgrid(t,y);
dT = ones(size(T));
dY = -3*Y;
figure('Name','impeule0r N=6')
quiver(T,Y,dT,dY)
axis tight
hold on
7
Lab 3 - Nathan Flo-
rean - MAT 275 Lab
Part (b)
t = linspace(0,5,100);
y = 2*exp(-3*t);
plot(t,y,'k-','linewidth',2);
8
Lab 3 - Nathan Flo-
rean - MAT 275 Lab
Part (c)
f=@(t,y)(-3*y);
clear t6 y6;
[t6,y6] = impeuler(f,[0,5],2,6);
plot(t6,y6,'ro-','linewidth',2);
hold off; % end plotting in this figure window
9
Lab 3 - Nathan Flo-
rean - MAT 275 Lab
The numerical solution when the stepsize N=6 is much more accurate when the improved Euler approxi-
mation is applied. When the improved Euler approximation is applied, there was no longer anything sig-
nificant when the stepsize N=6.
Part (d)
figure('Name','impeuler N=12');
% Plot slopefield for dy/dt = -3y.
t = 0:0.3:5; y = -1:0.4:2; % 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 = -3*Y; % dy = -3*y; this is the ODE
quiver(T,Y,dT,dY) % draw arrows (t,y)->(t+dt, t+dy)
axis tight % adjust look
hold on
t = linspace(0,5,100);
y = 2*exp(-3*t);
plot(t,y,'k-','linewidth',2);
t = 0:0.3:5;
y = -1:0.4:2;
f=@(t,y)(-3*y);
clear t12 y12;
[t12,y12] = impeuler(f,[0,5],2,12);
10
Lab 3 - Nathan Flo-
rean - MAT 275 Lab
plot(t12,y12,'ro-','linewidth',2);
hold off
Published with MATLAB® R2017b
11
