Question 1: Newton-Cotes quadrature

1

(a)    (i)        (x + 1)dx = 2

_1

1                                                                         344

(ii)    _1 ┌(2 - x)3 + (x - 2)2 ┐ dx =  12   (or 28.6666 . . .)

4                                                                         68062

(iii)         (7x6  - 3x4  - 9x2  - 5)dx = (or 13612.4)

3                                                                                 5

(b)  Modified quadrature functions:

function  q_mid=myMidpoint(f,a,b,N)

%  MYQUADALL  Numerical  qudrature  using  midpoint  rules

%

%  Usage :  q_mid=myMidpoint(f ,a ,b ,N)

%

%  Inputs :  f  =  function  handle  to  integrand

%                a ,b  =  limits  of  integration

%                N  =  number  of  panels  to  use

%

%  Outputs :  q_mid  =  midpoint  rule  approximation  to  the  integral

%  uniform  point  spacing

h=(b-a)/N;

%  datapoints

x=linspace(a+h/2,b-h/2,N);

%  integral

q_mid=h*sum(f(x));

end

function  q_simp=mySimpson(f,a,b,N)

%  MYQUADALL  Numerical  qudrature  using  Simpson✬s  rule

%

%  Usage :  q_simp=mySimpson(f ,a ,b ,N)

%

%  Inputs :  f  =  function  handle  to  integrand

%                a ,b  =  limits  of  integration

%                N  =  number  of  panels  to  use

%

%  Outputs :  q_simp  =  Simpson✬s  rule  approximation  to  the  integral

%  uniform  point  spacing

h=(b-a)/N;

%  datapoints

x=linspace(a,b,N+1);

%  integral  -  even  N  only

if  mod(N,2)==0

q_simp=(h/3)*(  f(a)  +  4*sum(f(x(2:2:end-1)))  +  2*sum(f(x(3:2:end-2)))  +  f(b)  );

else

q_simp=NaN;

end

end

Computed values of the integrals are as follows:

1

(i)        (x + 1)dx: _1

N	midpoint integral error		trapezium integral error		Simpson integral error
2〇 21 24 28	2 2 2 2	0 0 0 0	2 2 2 2	0 0 0 0	n/a 2 2 2	n/a 0 0 0

1

(ii) ┌(2 - x)3 + (x - 2)2 ┐ dx

_1

N	midpoint integral error		trapezium integral error		Simpson integral error
2〇 21 24 28	24.000000 27.500000 28.642857 28.666595	4.666667 1.166667 0.023810 0.000071	38.000000 31.000000 28.714286 28.666809	9.333333 2.333333 0.047619 0.000142	n/a 28.666667 28.666667 28.666667	n/a 0.000000 0.000000 0.000000

4

(iii)       (7x6 - 3x4 - 9x2 - 5)dx

3

N	midpoint integral error		trapezium integral error		Simpson integral error
2〇 21 24 28	12302.421875 13277.877686 13607.136903 13612.379439	1309.978125 334.522314 5.263097 0.020561	16264.500000 14283.460938 13622.926687 13612.441122	2652.100000 671.060938 10.526687 0.041122	n/a 13623.114583 13612.402628 13612.400000	n/a 10.714583 0.002628 0.000000

All three schemes have zero error for example (i), for all (allowed) N , while only Simpson’s rule has zero error for example (ii).  That is because the trapezium and midpoint rules have degree of accuracy 1, so they will exactly integrate polynomials of up to degree 1 (i.e. linear functions, like example (i)) with N = 1 (and, of course, all N > 1 too).  Simpson’s rule has degree of accuracy 3, so it will exactly integrate polynomials of up to degree 3 (i.e. cubics, like examples (i) and (ii)) with N = 2.

Note also the linear rate of convergence for the midpoint and trapezium rules, and the quadratic convergence of Simpson’s rule (where it is not exact). The midpoint rule also seems to have error

approximately half that of the trapezium rule (which is borne out by the error analysis).

(c)  Modified integration function:

function  q=myQuadIter(f,a,b,tol)

%  MYQUADITER  Numerical  qudrature  using  the  Trapezium  rule

%

%  Usage :  q=myQuadIter(f ,a ,b ,tol)

%

%  Inputs :  f  =  function  handle  to  integrand

%                a ,b  =  limits  of  integration

%                tol  =  desired  absolute  error  tolerance

%

%  Output :  q  =  approximation  to  the  integral

%  maximum  number  of  iterations

itmax  =  25;

%  set  intitial  values  of  iteration  count ,  subintervals ,  and  (dummy)  error

it  =  0;

N  =  2;

err  =  Inf ;

%  datapoints

x  =  linspace(a,b,N+1);

%  uniform  point  spacing

h  =  (b-a)/N;

%  first  approximation  to  the  integral

q  =  (h/2)*(f(a)+2*sum(f(x(2:end-1)))+f(b));

while  err>tol  &&  it

%  update  counters

q_old  =  q;

it  =  it+1;

N  =  2*N;

%  update  datapoints  and  current  integral  approx .

x  =  linspace(a,b,N+1);

h  =  (b-a)/N;

q  =  (h/2)*(f(a)+2*sum(f(x(2:end-1)))+f(b));

err  =  abs(q-q_old);

end

if  err>tol

error( ✬No  convergence  after  %d  steps ✬ ,it)

end

%  print  final  summary

fprintf(✬\nConvergence  after  %d  steps  (N=%d);  final  integral  value  is  %-22 .16g\n✬ , . . .

it,  N,  q);

fprintf( ✬Final  absolute  error  is  %g\n\n✬,  err);

end

The integral (iii) is around 105 , so to get 13 significant figures we need 8 decimal places of accuracy, and hence an absolute error tolerance of 5 × 10_9 . Calling the myQuadIter function

f  =  @(x)  7*x .^6  -  3*x .^4  -  9*x .^2  -  5;

q  =  myQuadIter(f,3,4,5e-9)

we get the integral correct to 13 s.f., with N = 2097152 (one very good reason why increasing N by 1 in the iteration is a bad idea!)

(d)  Modified integration function:

function  [qvec,evec,hvec]=myQuadAllIter(f,a,b,tol)

%  MYQUADALLITER  Numerical  qudrature  using  the  Trapezium  rule

%

%  Usage :  [qvec ,evec ,hvec]=myQuadAllIter(f ,a ,b ,tol)

%

%  Inputs :  f  =  function  handle  to  integrand

%                a ,b  =  limits  of  integration

%                tol  =  desired  absolute  error  tolerance

%

%  Output :  qvec  =  vector  of  approximations  to  the  integral

%

%

evec  =  vector  of  approximate  absolute  errors

hvec  =  vector  of  values  of  h

%  maximum  number  of  iterations

itmax  =  25;

%  set  intitial  values  of  iteration  count ,  subintervals ,  and  (dummy)  error

it  =  0;

N  =  2;

err  =  Inf ;

%  datapoints

x  =  linspace(a,b,N+1);

%  uniform  point  spacing

h  =  (b-a)/N;

%  first  approximation  to  the  integral

q  =  (h/2)*(f(a)+2*sum(f(x(2:end-1)))+f(b));

%  initialise  output  variables

qvec  =  q;

evec  =  err;

hvec  =  h;

while  err>tol  &&  it

%  update  counters

q_old  =  q;

it  =  it+1;

N  =  2*N;

%  update  datapoints  and  current  integral  approx .

x  =  linspace(a,b,N+1);

h  =  (b-a)/N;

q  =  (h/2)*(f(a)+2*sum(f(x(2:end-1)))+f(b));

err  =  abs(q-q_old);

%  update  output  variables

qvec  =  [qvec  q];

evec  =  [evec  err];

hvec  =  [hvec  h];

end

if  err>tol

error( ✬No  convergence  after  %d  steps ✬ ,it)

end

%  print  final  summary

fprintf(✬\nConvergence  after  %d  steps  (N=%d);  final  integral  value  is  %-22 .16g\n✬ , . . .

it,  N,  q);

fprintf( ✬Final  absolute  error  is  %g\n\n✬,  err);

end

Commands to plot a log-log graph of the absolute error versus h:

f  =  @(x)  7*x .^6  -  3*x .^4  -  9*x .^2  -  5;

[qvec,evec,hvec]  =  myQuadAllIter(f,3,4,5e-9);

loglog(hvec,evec, ✬r- ✬ , ✬LineWidth✬,2);

grid  on;

xlabel(✬h✬);  ylabel( ✬absolute  error✬);

xlim([1e-7  1e0]);  ylim([1e-9  1e3]);

100

10-5

10-6

10-4

h

10-2

100

The graph clearly shows a straight line, which tells us that error x hq, where q is the slope of the graph. Reading off the graph, the slope q = 2, and so the trapezium rule is a second order method.

Question 2: Numerical differentiation

(a)  The following commands will find numerical approximations to the derivative, using forward and

central difference, as well as the absolute errors:

h=2 .^(- (1:10));  x=1;

f=@(x)  cos(x);  fdash=@(x)  -sin(x);

fdashfd=(f(x+h)-f(x)) ./h;

fdashcd=(f(x+h)-f(x-h)) ./(2*h);