**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 m^{2}

and *Q*=450 m^{3}/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 m^{2}, *Q*=450 m^{3}/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);