b)
We use the linear transformation
x = (b-a)t+a.
I_a^b f(x) dx = I_0^1 f((b-a)t+a)(b-a)dt
~= (b-a)[Af(a) + Bf(1/4(b+3a)) + Cf(1/2(b+a)) + Df(1/4(3b+a)) + Ef(b)]

c)
ln(2) = I_0^1 1/(x+1) dx.
Using a) this is 0.6931746.
ln(2) = 0.6931471.
So our absolute error is 2.74x10^{-5}

2
a)
Let x = u^2. We have dx = 2u du.
I = I_0^4 1/(2x+8S(x)) dx, S=sqrt,
= I_0^2 2u/[2(u^2+4u)] du
= I_0^2 1/(u+4) du
= [ln(u+4)]_0^2 = ln(6) - ln(4) = ln(3/2)
= 0.40546510810816

b)
We want |f''(n)(1/12)h^2(b-a)| <= 10^-5.
We take f(x) = 1/(x+4) and a = 0, b = 2.
This avoids the singular evaluation at x = 0.
f''(x) = 2/(x+4)^3 and is decreasing in [0, 2].
So |f''(n)| <= 2/4^3.
Hence we need
2/4^3*1/12*h^2*2 <= 10^-5
-> h^2 <= 0.00192 -> h <= 0.0438 -> s = (b-a)/h >= 45.7~=46.

function Q = Trap(fname, s, a, b)
h = (b-a)/s;
Q = 0;
 for i= 1: s-1
    g = feval(fname, a+i*h);
    Q = Q + g;
end
 g1 = feval(fname, a);
 g2 = feval(fname, b);
 Q = h/2*(g1 +g2 + 2*Q);

function y = q2(x)
y = 1/(x+4);

>> Trap('q2',46,0,2)

ans =

0.40547057780408

3

a)
The error for Compound Simpson's rule with s panels, Q_s, is
E_s = I-Q_s ~= Ch^4, where C is some constant, h = (b-a)/s,
I = true integral value.
So, E_2s = I-Q_2s ~= Ch^4/16.
If we note that E_s ~= 16E_2s and subtract the 2 errors to cancel I, we have
E_2s ~= 1/15(Q_2s - Q_s).

function [Q, E] = SimpleSimp(fname, a, b)

h = (b-a)/2;
x1 = a+h/2;
x2 = a+h;
x3 = a+3*h/2;

fa = feval(fname, a);
f1 = feval(fname, x1);
f2 = feval(fname, x2);
f3 = feval(fname, x3);
fb = feval(fname, b);

S1 = h/3*(fa+4*f2+fb);
S2 = (h/6)*(fa + 4*f1 + 2*f2 + 4*f3 + fb);
Q = S2;
E = (1/15)*(S2-S1);

b)
function Q = AdaptSimp(fname,a,b,tol,lv)
% Implementation of the adaptive Simpsons rule.
% It uses at most 10 recursion levels.

[Q,E] = SimpleSimp(fname,a,b);
fprintf('Q = %f, E = %f, tol = %f, lv = %d, [%f, %f]\n', Q,E,tol,lv,a,b);
if((abs(E) > tol) & (lv <= 10))
    m = (a+b)/2;
    lv = lv+1;
    Q = AdaptSimp(fname,a,m,tol/2,lv) + AdaptSimp(fname,m,b,tol/2,lv);
end

function y = f(x)
y = (1-x)^.25;

>> AdaptSimp('f',0,1,.001,1)
Q = 0.769387, E = 0.002808, tol = 0.001000, lv = 1, [0.000000, 1.000000]
Q = 0.463639, E = 0.000002, tol = 0.000500, lv = 2, [0.000000, 0.500000]
Q = 0.323487, E = 0.001181, tol = 0.000500, lv = 2, [0.500000, 1.000000]
Q = 0.194936, E = 0.000001, tol = 0.000250, lv = 3, [0.500000, 0.750000]
Q = 0.136010, E = 0.000496, tol = 0.000250, lv = 3, [0.750000, 1.000000]
Q = 0.081961, E = 0.000000, tol = 0.000125, lv = 4, [0.750000, 0.875000]
Q = 0.057185, E = 0.000209, tol = 0.000125, lv = 4, [0.875000, 1.000000]
Q = 0.034460, E = 0.000000, tol = 0.000063, lv = 5, [0.875000, 0.937500]
Q = 0.024043, E = 0.000088, tol = 0.000063, lv = 5, [0.937500, 1.000000]
Q = 0.014489, E = 0.000000, tol = 0.000031, lv = 6, [0.937500, 0.968750]
Q = 0.010109, E = 0.000037, tol = 0.000031, lv = 6, [0.968750, 1.000000]
Q = 0.006092, E = 0.000000, tol = 0.000016, lv = 7, [0.968750, 0.984375]
Q = 0.004250, E = 0.000016, tol = 0.000016, lv = 7, [0.984375, 1.000000]

ans =

0.79982683519742
 

The plot of f(x) = (1-x)^1/4 shows that f(x) is smooth from 0 to 0.5 but
drops very sharply near 1. So we would expect that the integral will be
easy to approximate from 0 to 0.5 but difficult to approximate near 1.
From the many refinements near 1 and no refinement in [0, .5],
our output confirms our prediction.

b)
I_0^1 (1-x)^.25 dx
= 4/5[(1-x)^(5/4)]_0^1 = 4/5 = 0.8
Our absolute error is 0.000173 <= 0.001.
Computed solution is accurate within tolerance.

4
function Q = dbquad(fname, a, b, n)
%
% Computes the double integral I_a^b I_x^2x f(x,y) dy dx
% by using Compound Simpson's rule with n panels in both directions
% Here fname is the function name.
%

h = (b-a)/n;
w = 0;
u = 0;
for i= 1: n-1
    x1 = a+i*h;
    ay = g1(x1);
    by = g2(x1);
    u = u + phi(fname, x1, ay, by, n);
end
for i= 1: n
    x2 = a+(i-1/2)*h;
    ay = g1(x2);
    by = g2(x2);
    w = w + phi(fname, x2, ay, by, n);
end

ay = g1(a);
by = g2(a);
ga = phi(fname, a, ay, by, n);
ay = g1(b);
by = g2(b);
gb = phi(fname, b, ay, by, n);

Q = h/6*(ga + gb + 2*u +4*w );

function I = phi(fname, x, a, b, n)
%
% Compound Simpson's rule with n panels in the y-direction.
%
h = (b-a)/n;
w = 0;
u = 0;
for i= 1: n-1
    u = u + feval(fname, x, a+i*h);
    w = w + feval(fname, x, a+(i-1/2)*h);
end
ga = feval(fname, x, a);
gb = feval(fname, x, b);
g = feval(fname, x, a+(n-1/2)*h);
I = h/6*(ga + gb + 2*u +4*(w + g));

function y = g1(x)
y = x;

function y = g2(x)
y = 2*x;

function z = q4func(x,y)
z = y^2+x^3;
 

>> Q = dbquad('q4func', 0, 1, 6)

Q =

0.78333976337449

I_0^1 I_x^2x x^3+y^2 dy dx
= I_0^1 x^4+7/3x^3 dx
= 47/60
= 0.78333333333333
Our absolute error is
6.43e-06.

Bonus
The orthonormal set is given by:
{q0, q1, q2, q3} = {1/S(2), S(3/2)x, S(45/8)(x^2-1/3), S(175/8)(x^3-3/5x)}
where S denotes the sqrt function.

Let f(x) = cos(pi x)
p*(x) = <f, q0>q0 + <f, q1>q1 + <f, q2>q2 + <f, q3>q3

<f, q0>
= I_{-1}^1 1/S(2) cos(pi x) dx
= 1/(S(2)pi)[sin(pi x)]_{-1}^1 = 0

<f, q1> and <f, q3> = 0 since their integrands are odd functions in [-1, 1].

<f, q2>
= S(45/8)(-4/pi^2) by using integration by parts.
So p*(x) = (-4/pi^2)(45/8)(x^2-1/3).
Plot