MAT275 MATLAB 3 Buyun Ma
Exercise 1
Part(a)
f = @(t,y)(1.5*y);
t = linspace(0,1.4,100);y=-2*exp(1.5*t);
[t7,y7]=euler(f,[0,1.4],-2,7);
[t70,y70]=euler(f,[0,1.4],-2,70);
[t700,y700]=euler(f,[0,1.4],-2,700);
[t7000,y7000]=euler(f,[0,1.4],-2,7000);
y7(end);
y70(end);
y700(end);
y7000(end);
e7 = y(end)-y7(end);
e70 = y(end)-y70(end);
e700 = y(end)-y700(end);
e7000 = y(end)-y7000(end);
r7 = e7/e70;
r70 = e70/e700;
r700 = e700/e7000;
>> disp('|N |approximation| error | ratio |')
disp('|7 |-12.5497 |-3.7826| N/A |')
disp('|70 |-15.8356 |-0.4967|7.6156 |')
disp('|700 |-16.2811 |-0.0513|9.6891 |')
disp('|7000|-16.3272 |-0.0051|9.9679 |')
output:
_______________________________________
|N |approximation| error | ratio |
|7 |-12.5497 |-3.7826| N/A |
|70 |-15.8356 |-0.4967|7.6156 |
|700 |-16.2811 |-0.0513|9.6891 |
|7000|-16.3272 |-0.0051|9.9679 |
_______________________________________
Part(b):
As ‘N’ increases the ratio of consecutive errors decreases. This states that as the number of steps
are increased, Euler method tries to minimize the error (as the ratio of consecutive errors
increases). Also a quick look at column 3 (Error) suggests that as the step size is decreased from
7000 to 700, error also is also changed by the same factor (here 7000/700). This indicates that
the Euler method is of order ‘h’ that is everytime the step size is decreased by a factor(k), error
also changes by the same factor (k).
Part(c):
Euler's method underestimates the actual curve, because the concavity of the function is facing
upwards, so Euler's will always under estimate.
Exercise 2
(a)
t = 0:0.2:3; y = -8:1.5:7;
[T,Y] = meshgrid(t,y);
dT = ones(size(T));
dY = -5.1*Y;
quiver(T,Y,dT,dY)
axis tight
hold on
(b)
t = linspace(0,3,100);
y = 2*exp(-5.1*t);
plot(t,y,'k','LineWidth',2);
(c)
f = inline('-5.1*y','t','y');
[T,Y] = euler(f,[0,3],2,7);
plot(T,Y,'ro-','Linewidth',2);
(d)
figure
t = 0:0.2:3;
y = -1.2:0.4:2.6;
[T,Y] = meshgrid(t,y);
dT = ones(size(T));
dY = -5.1*Y;
quiver(T,Y,dT,dY);
axis tight;
hold on;
t = linspace(0,3,100);
y = 2*exp(-5.1*t);
plot(t,y,'k','LineWidth',2);
f = inline('-5.1*y','t','y');
[T,Y] = euler(f,[0,3],2,14);
plot(T,Y,'ro-', 'LineWidth',2);
Exercise 3
(impeuler.m)
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
m1=f(t(n),y(:,n));
m2=f(t(n+1),y(:,n)+h*m1);
y(:,n+1)=y(:,n)+h*(m1+m2)/2;
end
t=t';
y=y';
end
f = inline('1.5*y','t','y');
[t7,y7] = impeuler(f,[0,1.4],-2,7);
[t7,y7]
ans =
0 -2.0000
0.2000 -2.6900
0.4000 -3.6180
0.6000 -4.8663
0.8000 -6.5451
1.0000 -8.8032
1.2000 -11.8403
1.4000 -15.9252
Exercise 4
f = inline('1.5*y','t','y');
t = linspace(0,1.4,100);
y = -2*exp(1.5*t);
[t7,y7] = impeuler(f,[0,1.4],-2,7);
[t70,y70] = impeuler(f,[0,1.4],-2,70);
[t700,y700] = impeuler(f,[0,1.4],-2,700);
[t7000,y7000] = impeuler(f,[0,1.4],-2,7000);
y7(end);
y70(end);
y700(end);
y7000(end);
e7 = y(end) - y7(end);
e70 = y(end) - y70(end);
e700 = y(end) - y700(end);
e7000 = y(end) - y7000(end);
r7 = e7/e70;
r70 = e70/e700;
r700 = e700/e7000;
output:
N
approximation
error
ratio
7
-15.9252
-0.4071
N/A
70
-16.3273
-0.0050
80.9417
700
-16.3323
-5.1331e-05
97.9823
7000
-16.3323
-5.1435e-07
99.7976
As ‘N’ increases the ratio of consecutive errors also increases. This states that as the number of steps
are increased, Improved Euler method tries to minimize the error (as the ratio of consecutive errors
increases). Also a quick look at column 3 (Error) suggests that as the step size is increased from 7000
to 700, error is approximately changed by the square of the factor (here 7000/700). This indicates that
the Improved Euler method is of order ‘h2’ that is everytime the step size is decreased by a factor(k),
error changes by the square of the factor (k).
Exercise 5
(a)
figure
t = 0:0.2:3;
y = -1.2:0.4:2.6;
[T,Y] = meshgrid(t,y);
dT = ones(size(T));
dY = -5.1*Y;
quiver(T,Y,dT,dY);
axis tight;
hold on;
(b)
t = linspace(0,3,100);
y = 2*exp(-5.1*t);
plot(t,y,'k','LineWidth',2);
(c)
f = inline('-5.1*y','t','y');
[T,Y] = impeuler(f,[0,3],2,7);
plot(T,Y,'ro-','Linewidth',2);
(d)
figure
t = 0:0.2:3;
y = -1.2:0.4:2.6;
[T,Y] = meshgrid(t,y);
