Autumn Schuelka Wednesday 8:35am
Autumn Schuelka
Wednesday 8:35am
Palou-Genevive Toutain
MAT 275 MATLAB Lab 3
Exercise 1
% This is the euler.m function
function [t,y] = euler(f,tspan,y0,N) % Solves the IVP y’=f(t,y,y(t0)=y0
over the time interval
tspan=[t0,tf]
m = length(y0);
t0 = tspan(1); % t0=initial time value
tf = tspan(2); % tf=final time value
h = (tf-t0)/N; % N=number of steps used, h evaluates step size
t = linspace(t0,tf,N+1);
y = zeros(m,N+1);
y(:,1) = y0'; % y0=initial value of dependent variable
for n = 1:N
y(:,n+1) = y(:,n)+h*f(t(n),y(:,n)); % implements Euler’s method
end
t = t'; y = y';
end
(a)
f = inline('2*y','t','y') % defines the function y’=2y
f =
Inline function:
f(t,y) = 2*y
t = linspace(0,.5,100); y = 3*exp(2*t); % evaluates the exact solution
[t50,y50] = euler(f,[0,.5],3,50); % evaluates when N=50
[t50,y50]
ans =
0 3.0000
: :
0.5000 8.0748 % edited the output to only show first and last entries
[t500,y500] = euler(f,[0,.5],3,500); % solves the ODE using Euler when N=500
[t,y] = euler(f,[0,.5],3,500);
[t,y]
ans =
0 3.0000
: :
0.5000 8.1467
[t5000,y5000] = euler(f,[0,.5],3,5000); % solves with 5000 steps
[t,y] = euler(f,[0,.5],3,5000);
[t,y]
ans =
0 3.0000
: :
0.5000 8.1540
Autumn Schuelka Wednesday 8:35am
e500 = y(end)-y500(end) % error with N = 500
e500 =
0.0081
e5000 = y(end)-y5000(end) % error with N = 5000
e5000 =
8.1534e-04
ratio = e50/e500 % ratio of error between N = 50 and N = 500
ratio =
9.8381
ratio = e500/e5000 % ratio of error between N = 500 and N = 5000
ratio =
9.9835
N
Approximation
Error
Ratio
5
7.4650
0.6899
--
50
8.0748
0.0801
8.6148
500
8.1467
0.0081
9.8381
5000
8.1540
8.1534e-4
9.9835
(b)
As the number of steps used increases, the ratio of consecutive errors
increases. As the number of steps increases by a factor of 10, the error
decreases by a factor of 10. There is a theorem that states if y’=f(t,y),
and y(t0)=y0 has a unique solution y(t) on the closed interval [t0,tk] and
y(t) has a continuous second order derivative on this interval. Then there
exists a constant C such that the following is true:
the approximations compared to the actual values computed using Euler’s
method where stepsize h>0, then
|y(tn)-yn| Ch for each n = 1,2,3,…,k
(c)
From the geometrical representation of Euler’s method, the tangent line is
used to determine the next value via the derivative. Since the slope of the
actual value graph is constantly changing, the tangent line is only a single
point of the whole graph over each interval, or stepsize. Using these values
makes the Euler approximation underestimate because the slope of the actual
value graph is changing at a greater rate than the tangent line.
Autumn Schuelka Wednesday 8:35am
Exercise 2
(a)
t = 0:.45:10; y = -30:6:42; % define grid of values in t and y direction
[T,Y] = meshgrid(t,y); % creates d 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) % draws arrows (t,y)->(t+dt, t+dy)
axis tight % adjust look
hold on
% output graph of vfield
(b)
t = linspace(0,10,200);
y = 3*exp(2*t);
plot(t,y,'k-','linewidth',2)
% output graph
Autumn Schuelka Wednesday 8:35am
(c)
f = inline('-2*y','t','y')
f =
Inline function:
f(t,y) = -2*y
[t,y] = euler(f,[0,10],3,8);
plot(t,y,'ro-','linewidth',2)
% output graph that is just like Figure L3b
The approximations for the values are so inaccurate because they are not
determined from the original function, but the slope of the tangent line of
that function at the particular points.
Autumn Schuelka Wednesday 8:35am
(d)
figure
t = 0:.4:10; y = -1:0.4:3;
[T,Y] = meshgrid(t,y);
dT = ones(size(T));
dY = -2*Y;
quiver(T,Y,dT,dY)
axis tight
hold on
t = linspace(0,10,200); y = 3*exp(-2*t);
plot(t,y,'k-','linewidth',2)
[Insert graph]
f = inline('-2*y','t','y')
f =
Inline function:
f(t,y) = -2*y
[t,y] = euler(f,[0,10],3,16);
plot(t,y,'ro-','linewidth',2)
% output graph that is just like Figure L3c
Autumn Schuelka Wednesday 8:35am
Exercise 3
% This is the impeuler.m function
function [t,y] = impeuler(f,tspan,y0,N)
m = length(y0);
t0 = tspan(1);
tf = tspan(2);
h = (tf-t0)/N;
t = linspace(t0,tf,N+1);
y = zeros(m,N+1);
y(:,1) = y0';
for n = 1:N
f1 = f(t(n),y(:,n));
f2 = f(t(n+1), y(:,n)+h*f1);
y(:,n+1) = y(:,n)+h/2*(f1+f2);
end
t = t'; y = y';
end
f = inline('2*y','t','y') % defines the ODE y’=2y
f =
Inline function:
f(t,y) = 2*y
[t5,y5] = impeuler(f,[0,.5],3,5); % computes the solution when N=5
[t5,y5]
ans =
0 3.0000
0.1000 3.6600
0.2000 4.4652
0.3000 5.4475
0.4000 6.6460
0.5000 8.1081
The output values are the same as in the Lab instructions.
Autumn Schuelka Wednesday 8:35am
Exercise 4
(a)
[t50,y50] = impeuler(f,[0,.5],3,50); % computes when N=50
[t50,y50]
ans =
0 3.0000
: :
0.5000 8.1543 % edited the output to only show first and last entries
[t500,y500] = impeuler(f,[0,.5],3,500); % computes when N=500
[t500,y500]
ans =
0 3.0000
: :
0.5000 8.1548
[t5000,y5000] = impeuler(f,[0,.5],3,5000); % computes when N=5000
[t5000,y5000]
ans =
0 3.0000
: :
0.5000 8.1548
t = linspace(0,.5,100); y = 3*exp(2*t); % defines exact solution of the ODE
e5 = y(end)-y5(end) % error with N = 5 % computes error N=5
e5 =
0.0467
e50 = y(end)-y50(end) % error with N = 50 % computes error N=50
e50 =
5.3555e-04
e500 = y(end)-y500(end) % error with N = 500 % computes error N=500
e500 =
5.4284e-06
e5000 = y(end)-y5000(end) % error with N = 5000 % computes error N=5000
e5000 =
5.4357e-08
ratio = e5/e50 % ratio of error between N = 5 and N = 50
% computes ratio between N=5 and N=50
ratio =
87.2394
ratio = e50/e500 % ratio of error between N = 50 and N = 500
% computes ratio between N=50 and N=500
ratio =
98.6567
ratio = e500/e5000 % ratio of error between N = 500 and N = 5000
% computes ratio between N=500 and N=5000
ratio =
99.8650
N
Approximation
Error
Ratio
5
8.1081
0.0467
--
50
8.1543
5.3555e-4
87.2394
500
8.1548
5.4284e-6
98.6567
5000
8.1548
5.4357e-8
99.8650
Autumn Schuelka Wednesday 8:35am
(b)
Very similarly to the example in exercise 1, As the number of steps used
increases, the ratio of consecutive errors increases. As the number of steps
increases, the error decreases. There is a theorem that states if y’=f(t,y),
and y(t0)=y0 has a unique solution y(t) on the closed interval [t0,tk] and
y(t) has a continuous third order derivative on this interval. Then there
exists a constant C such that the following is true:
the approximations compared to the actual values computed using Improved
Euler’s method where stepsize h>0, then
|y(tn)-yn| Ch2 for each n = 1,2,3,…,k
Exercise 5
(a)
t = 0:.45:10; y = -30:6:42; % define grid of values in t and y direction
[T,Y] = meshgrid(t,y); % creates d 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) % draws arrows (t,y)->(t+dt, t+dy)
axis tight % adjust look
hold on
Autumn Schuelka Wednesday 8:35am
(b)
t = linspace(0,10,200);
y = 3*exp(2*t);
plot(t,y,'k-','linewidth',2)
(c)
f = inline('-2*y','t','y')
f =
Inline function:
f(t,y) = -2*y
[t,y] = impeuler(f,[0,10],3,8);
plot(t,y,'ro-','linewidth',2)
Autumn Schuelka Wednesday 8:35am
(d)
figure
t = 0:.4:10; y = -1:0.4:3;
[T,Y] = meshgrid(t,y);
dT = ones(size(T));
dY = -2*Y;
quiver(T,Y,dT,dY)
axis tight
hold on
t = linspace(0,10,200); y = 3*exp(-2*t);
plot(t,y,'k-','linewidth',2)
[Insert graph]
f = inline('-2*y','t','y')
f =
Inline function:
f(t,y) = -2*y
[t,y] = impeuler(f,[0,10],3,16);
plot(t,y,'ro-','linewidth',2)
The results of Improved Euler’s method compared to the original Euler’s
method show that Improved Euler’s method is more accurate at smaller step
size intervals, when you begin to increase the step size the approximation
is less accurate, because h is squared. Overall, Euler’s method is not
entirely accurate, but when the step size is increased, the error is only
changed by the inverse of the factor of step size. This makes it more
accurate at larger step size values.
