Class Code:

%MATLAB

% Romberg and Gauss Integration methods

% function to integrate over [0,2]

f = @(x) exp(x).*sin(x)./(1+x.^2);

% make a plot first

x = 0:0.01:2;

y = f(x);

plot(x,y);

% trapezoid rule estimate (what h needed to converge?)

% (there is no analytical solution...)

a = 0; b = 2; n = 2000;

I = trapezoid(f,a,b,n);

disp(['Trapezoid rule: ', num2str(n), ' fxn evals, estimate: ' num2str(I,

16)]);

% Built-in Matlab: Adaptive quadrature

[I, num_evals] = quad(f,a,b);

disp(['Built-in "quad": ', num2str(num_evals), ' fxn evals, estimate: '

num2str(I, 16)]);

% Romberg estimate

num_rows = 4;

R = romberg(f,a,b,num_rows);

disp(R) % show the whole table

I = R(1,num_rows+1); % the best integral estimate

num_evals = sum(2.^(1:num_rows+1));

disp(['Romberg: ', num2str(num_evals), ' fxn evals, estimate: ' num2str(I,

16)]);

% Gauss Quadrature estimates

disp(['Gauss 2-pt estimate: ' num2str(gauss2(f,a,b), 16)]);

disp(['Gauss 3-pt estimate: ' num2str(gauss3(f,a,b), 16)]);

disp(['Gauss 4-pt estimate: ' num2str(gauss4(f,a,b), 16)]);

% Functions used above: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

% Trapezoid rule: function, a,b, number of segments

function I = trapezoid(f,a,b,n)

h = (b-a);

I = (h/2)*(f(a) + 2*sum(f(a+h:h:b-h)) + f(b));

end

% Romberg: function, a, b, number of rows in the table

% Follows Figure 22.4 of Chapra & Canale 7th Ed.

% Returns the full Romberg table, R

function R = romberg(f,a,b,num_rows)

R = zeros(num_rows, num_rows); % Romberg table

n = 1; % first step size: (b-a), one segment

R(1,1) = trapezoid(f,a,b,n); % first entry in table

% add a new row every iteration, and update columns

for iter = 1:num_rows

n = n*2; % every new row uses 2x the number of segments

R(iter+1,1) = trapezoid(f,a,b,n); % new row entry, bottom of first

column

% now that we've added a new row in column 1, update all other

% columns (one new value for each column)

for k = 2:iter+1

j = 2 + iter - k; % indexing hack to place the new value

R(j,k) = (4^(k-1) * R(j+1,k-1) - R(j,k-1)) / (4^(k-1) - 1);

end

end

end

% Gauss Quadrature (different functions for # of points)

% weights and x values from Table 22.1 of Chapra & Canale 7th

function I = gauss2(f,a,b)

c = [1 1]; % weight coefficients

x = [-1/sqrt(3) 1/sqrt(3)];

y = (a+b)/2 + (b-a)/2 * x;

I = (b-a)/2 * sum(c.*f(y));

end

function I = gauss3(f,a,b)

c = [0.5555556 0.8888889 0.5555556]; % weight coefficients

x = [-0.774596669 0 0.774596669];

y = (a+b)/2 + (b-a)/2 * x;

I = (b-a)/2 * sum(c.*f(y));

end

function I = gauss4(f,a,b)

c = [0.3478548 0.6521452 0.6521452 0.3478548]; % weight coefficients

x = [-0.861136312 -0.339981044 0.339981044 0.861136312];

y = (a+b)/2 + (b-a)/2 * x;

I = (b-a)/2 * sum(c.*f(y));

end

22.15 The amount of mass transported via a pipe over a period of time can be computed as where M = mass(mg), ti = the initial time(min), t2 = the final time (min), 2(t) - flow rate (m*/min), and c(t) -concentration (mg/m*). The following functional representations define the tem poral variations in flow and concentration: Q(t) = 9 + 5 cos(0.41) c(t) = 5e-0.5, + 2e0.15, Determine the mass transported between2 and t2- 8 min with Romberg integration to a tolerance of 0.1% Repeat 22.15 using 5-point and 6-point Gauss quadrature. Compare to your answer from Problem 1, and also compare to the result using the built-in quad function