EMAT20920: Numerical Methods in MATLAB WORKSHEET 5 — SOLUTIONS

联系我们： 手动添加方式: 微信>添加朋友>企业微信联系人>13262280223 或者 QQ: 1483266981

EMAT20920: Numerical Methods in MATLAB

WORKSHEET 5 — SOLUTIONS

NUMERICAL CALCULUS

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) (7×6 – 3×4 – 9×2 – 5)dx = (or 13612.4)

3 5

(b) Modi ed 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) (7×6 – 3×4 – 9×2 – 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) Modi ed 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 .16gn , . . .

it, N, q);

fprintf( Final absolute error is %gnn , err);

end

The integral (iii) is around 105 , so to get 13 signi cant gures 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) Modi ed 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 .16gn , . . .

it, N, q);

fprintf( Final absolute error is %gnn , 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 o the graph, the slope q = 2, and so the trapezium rule is a second order method.

Question 2: Numerical di erentiation

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

central di erence, 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);

abserrfd=abs(fdashfd-fdash(x));

abserrcd=abs(fdashcd-fdash(x));