Euler Heun

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.

  1. Analytically, plotted with increments of 0.01
  2. Euler’s method with a step size of 0.5
  3. Midpoint method with a step size of 0.5
  4. 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)

  1. 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.

  1. 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);