Euler Heun
TASK 1
Consider the following ordinary differential equation (ODE):
Find y(3) using pen and paper with a step size of h=1, and the initial condition of y(0)=3 for the following methods:
A. Euler’s method
B. Midpoint method
C. Heun’s method
TASK 2
Write function files for the following methods:
A. Euler’s method
B. Midpoint method
C. Heun’s method
The function headers should look similar to the following:
function [t, y] = euler(dydt, tspan, y0, h)
function [t, y] = heun(dydt, tspan, y0, h)
function [t, y] = midpoint(dydt, tspan, y0, h)
TASK 3
The amount of a uniformly distributed radioactive contaminant contained in a closed reactor is
measured by its concentration c (Becquerel/litre or Bq/L). The contaminant decreases at a decay rate
proportional to its concentration; that is decay rate = –kc, where k is a constant with units of day-1. Solve
for the concentrations from t=0 to t=5 days with k=0.175 day-1 and a step size of 0.1 days using:
A. Euler’s method
B. Heun’s method
C. Midpoint method
D. MATLAB’s ode45
Assume the concentration at t=0 is 100 Bq/L. Plot the solutions of all parts in the one figure and print
out the concentration value at t=5 days for all parts.
TASK 4
Solve the following equation over the interval from t=0 to 2 where y(0) = 1 and display your results on
the same graph using the methods below.
- Analytically, plotted with increments of 0.01
- Euler’s method with a step size of 0.5
- Midpoint method with a step size of 0.5
- Ode45 using adaptive step size
Run your code multiple times with various step sizes. As you decrease the step size, which method (Euler or midpoint) converges faster? Why is this? Remember to include a legend for your plot.
TASK 5
A storage tank (shown below) contains a liquid at depth y where y=0 when the tank is half full. Liquid is
withdrawn at a constant flow rate Q to meet demands. The contents are resupplied at a sinusoidal rate
3Qsin2(t)
- The change in volume can be written as
where A is the surface area and is assumed to be constant. Use Euler’s method to solve for the
depth y from t=0 to t=10 days with a step size of 0.1 days. The parameter values are A=1250 m2
and Q=450 m3/day. Assume that the initial condition is y=0. Plot the solution.
- For the same storage tank, suppose that the outflow is not constant but rather depends on the
For this case, the differential equation for the change in volume can be written as
Use Euler’s method to solve for the depth y from t=0 to 10 days with a step size of 0.1 days. The
parameter values are A=1250 m2, Q=450 m3/day and α=150. Assume that the initial condition is y=0. Plot this solution on the previous figure. Remember to include a legend and place it in the north-west corner of the figure.
Solution
comsimp.m
function I=comsimp(func,a,b,n) %function to evaluate composite
simpson’s integration
h=(b-a)/(n-1); %step size
x=a:h:b; %array of equally sized points
tp=0; %temp variable to store the summation
for j=2:2:size(x,2)-2
tp=tp+(4*func(x(j)))+(2*func(x(j+1)));
end
I=(h/3)*(tp+func(x(1))+func(x(2)));
end
comtrap.m
function I=comtrap(func,a,b,n) %function to evaluate composite trapezoidal integration
h=(b-a)/(n-1); %step size
x=a:h:b; %array of equally sized points
tp=0; %temp variable to store the trapezoidal summation
for j=1:1:size(x,2)-1
tp=tp+func(x(j))+func(x(j+1));
end
I=(h/2)*tp;
end
comtrap2.m
function I=comtrap2(x,y) %This function uses data sets instead of the function
tp=0; %temporary variable for summation
for j=1:1:size(y,2)-1
tp=tp+(0.5*(x(j+1)-x(j))*(y(j)+y(j+1)));
end
I=tp;
end
ft2.m
function y=ft2(x) %function file for task 2
y=1+exp(-x);
end
ft3.m
function y=ft3(x) %function file for task 3
y=5+(0.25.*x.*x);
end
ft4.m
function y=ft4(x) %function file for task 4
y=200*(x/(5+x))*exp((-2*x)/30);
end
lab9t2.m
%Script to perform task 2
t2a=comtrap(@ft2,0,4,7); %implementing composite trapezoidal rule with 7 points
t2b=comsimp(@ft2,0,4,7); %implementing composite simpson’s rule with 7 points
fileID=fopen(‘lab9t2.txt’,’w’);
fprintf(fileID,’The value of the integral obtained by composite trapezoidal method using 7 points is %4.6f\r\n’,t2a);
fprintf(fileID,’The value of the integral obtained by composite simpsons method using 7 points is %4.6f\r\n’,t2b);
fclose(fileID);
lab9t3.m
%Script for task 3
t3a=integral(@ft3,0,11); %Analytical Integral
t3b=comtrap(@ft3,0,11,13); %Composite Trapezoidal Integral with n=13
t3c=comsimp(@ft3,0,11,13); %Composite Simpson’s Integral with n=13
fileID=fopen(‘lab9t3.txt’,’w’);
fprintf(fileID,’The value of the integral obtained analytically is %4.6f\r\n’,t3a);
fprintf(fileID,’The value of the integral obtained by composite trapezoidal method using 13 points is %4.6f\r\n’,t3b);
fprintf(fileID,’The value of the integral obtained by composite simpsons method using 13 points is %4.6f\r\n’,t3c);
fclose(fileID);
lab9t4.m
%Script for task 4
z=0:30/100:30; %equally spaced points between z=0 and z=H=30
fz=ft4(z); %values of the function f(z)
figure(1)
plot(z,fz,’rx-‘),xlabel(‘elevation above the deck’),ylabel(‘force on the sailboat’),title(‘Force vs. Elevation’),grid on
saveas(gcf,’lab9t4.png’,’png’); %saves plot to file
t4b=integral(@ft4,0,30); %analytical integral
n1=3;
n2=3;
t4c1=comtrap(@ft4,0,30,n1); %initial trapezoidal integral with n=3
t4c2=comsimp(@ft4,0,30,n2); %initial simpson integral with n=3
while abs((t4c1-t4b)/t4b)>0.01 %maximum error criterion
n1=n1+1; %decreasing step size with each check
t4c1=comtrap(@ft4,0,30,n1);
end
while abs((t4c2-t4b)/t4b)>0.01 %maximum error criterion
n2=n2+1; %decreasing step size with each check
t4c1=comtrap(@ft4,0,30,n2);
end
fileID=fopen(‘lab9t4.txt’,’w’);
fprintf(fileID,’The minimum number of points required to achieve 0.01% error using trapezoidal integral is %d\r\n’,t4c1);
fprintf(fileID,’The minimum number of points required to achieve 0.01% error using simpson integral is %d\r\n’,t4c2);
fclose(fileID);
lab9t5.m
%Script for task 5
y=[0 2 4 5 6 9];
H=[0.5 1.3 1.25 1.8 1 0.25];
U=[0.03 0.06 0.05 0.13 0.11 0.02];
Y=H.*U; %element-wise multiplication
figure(1)
plot(y,H,’rx-‘,y,Y,’ko-‘),xlabel(‘distance from the bank’),title(‘Data Plots’),legend(‘channel depth(H)’,’flow flux’),grid on
saveas(gcf,’lab9t5_1.png’,’png’);
t5a=polyfit(y,H,2); %second-order polynomial fit for H
t5b=polyfit(y,Y,2); %second-order polynomial fit for UH
t51=polyval(t5a,y); %evaluating polynomial fit
t52=polyval(t5b,y);
figure(2)
plot(y,H,’bo’,y,t51,’r-‘,y,Y,’bo’,y,t52,’k-‘),xlabel(‘distance from the bank’), title(‘Plots of polynomial fits for H(y) and U(y)H(y)’),grid on
saveas(gcf,’lab9t5_2.png’,’png’);
[email protected](x) t5a(1)*x*x+t5a(2)*x+t5a(3); %constructing polynomial splines
[email protected](x) t5b(1)*x*x+t5b(2)*x+t5b(3);
r5a=comtrap(ft51,0,9,7); %Composite Trapezoidal Integral on Polynomial Fit H with n=7
r5b=comtrap(ft52,0,9,7); %Composite Trapezoidal Integral on Polynomial Fit UH with n=7
r5c=comtrap2(y,H); %Integration by Trapezoidal Segments for H(y)
r5d=comtrap2(y,Y); %Integration by Trapezoidal Segments for U(y)H(y)
fileID=fopen(‘lab9t5.txt’,’w’);
fprintf(fileID,’The value of Ac obtained by composite trapezoidal method using 7 points is %4.6f\r\n’,r5a);
fprintf(fileID,’The value of Q obtained by composite trapezoidal method using 7 points is %4.6f\r\n’,r5b);
fprintf(fileID,’The value of Ac obtained by integrating trapezoidal segments is %4.6f\r\n’,r5c);
fprintf(fileID,’The value of Q obtained by integrating trapezoidal segments is %4.6f\r\n’,r5d);
fclose(fileID);