1 / 16100%
MAT 275 MATLAB LAB 3
Alejandro Carranza
Monday 9:40 AM
Javier Baez
Exercise 1
Part (a)
% Define ODE function f
f = inline('2*y','t','y');
% analytical solution vector.
t = linspace(0,.5,100); y=3*exp(2*t); % exact
solution
% Create vector
[t50,y50] = euler(f,[0,.5],3, 50); % solves the
ODE
% Solve IVP numerically with 5 time steps.
[t5,y5]=euler(f,[0,.5],3, 5);
% Solve IVP numerically 50 time steps.
[t50,y50]=euler(f,[0,.5],3, 50);
% Solve IVP numerically with 500 time steps.
[t500,y500]=euler(f,[0,.5],3, 500);
% Solve IVP numerically with 5000 time steps.
[t5000,y5000]=euler(f,[0,.5],3, 5000);
% Compute numerical with 5
% time steps.
e5 = y(end)-y5(end)
% define error
% Compute numerical solution error
e50 = y(end) - y50(end)
% Compute numerical solution error with 500
% time steps.
e500 = y(end)- y500(end)
% Compute numerical solution error with 5000
% time steps.
e5000 = y(end)- y5000(end)
% Compute ratio of errors
ratio = e5/e50
% Compute ratio of errors
ratio = e50/e500
% Compute ratio of errors between N=500 and N=5000.
ratio = e500/e5000
% Display table of errors
disp('----------------------------------------------')
disp('| N | approximation | error | ratio |')
disp('|------|-----------------|----------|--------|')
disp('| 5 | 7.4650 | 0.6899 | N/A |')
disp('| 50 |
disp('| 500 |
disp('| 5000 |
disp('----------------------------------------------')
8.0748 | 0.0801 | 8.6148 |')
8.1467 | 0.0081| 9.8381|')
8.154 | 8.1534e-04| 9.9835|')
e5 = 0.6899
e50 = 0.0801
e500 = 0.0081
e5000 =
8.1534e-04
ratio =
8.6148
ratio =
9.8381
ratio =
9.9835
----------------------------------------------
| N | approximation | error | ratio |
|------|-----------------|----------|--------|
| 5 | 7.4650 | 0.6899 | N/A |
| 50 |
| 500 |
| 5000 |
----------------------------------------------
8.0748 | 0.0801 | 8.6148 |
8.1467 | 0.0081| 9.8381|
8.154 | 8.1534e-04| 9.9835|
Part (b).As the step size is decreased the error gets smaller by the factor relative to the
step size taken. In the table you can see this as we change the step size from 5 to 50 the
error is reduced by a tenth.
Part (c). The reason that eulers method underestimates the function is because it finds
values relative to the slope of the function rather than the function itself.
Exercise 2
Part (a)
% Plot slope field for dy/dt = -2y.
t = 0:.45:10; y = -30:6:42 ; % define grid
[T,Y]=meshgrid(t,y); % creates 2d matrices
dT = ones(size(T)); % dt=1 for all points
dY = -2*Y; % dy = -2*y
quiver(T,Y,dT,dY) % draw arrows (t,y)->(t+dt, t+dy)
axis tight % adjust look
hold on
!
Part (b)
% Define vector t
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 slope field from
(a).
plot(t, y,'k', 'linewidth', 2)
Part (c)
% Define ODE function.
inline( '-2*y', 't', 'y')
% Compute numerical solution to IVP with 8
[t8,y8]=euler(f,[0,.5],3, 8);
% Plot numerical solution with analytical solution
plot(t8, y8, 'ro-','linewidth',2)
hold off; % end plotting in this figure window
ans =
!
Inline function:
ans(t,y) = -2*y
Part (d)
% Define new grid of t and y values
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
dY = -2*Y; % dy = -2*y
quiver(T,Y,dT,dY) % draw arrows (t,y)->(t+dt, t+dy)
axis tight % adjust look
hold on
% Plot slope field for dy/dt = -2y
% Define vector t
linspace(0,10, 100);
% Create vector of analytical solution
y = 3*exp(-2*y);
% Define ODE function.
inline('-2*y', 't', 'y')
% Compute numerical solution to IVP with 16 time
steps
[t16,y16]=euler(f,[0,.5],3, 16);
plot(t16,y16,'ro-','linewidth',2)
ans =
Inline function:
ans(t,y) = -2*y
Exercise 3
% Display contents
type 'impeuler.m'
% Define ODE function for dy/dt = f(t,y) = 2y.
inline('2*y', 't', 'y')
% Compute numerical solution to IVP with 5 time
steps