﻿ 数值分析作业九 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 数值分析作业九

Chap 5.1

Exercises

2. Use the three-point centered-difference formula to approximate f (0), where f(x) = e x , for (a) h = 0.1 (b) h = 0.01 (c) h = 0.001.

Solution:

Because: f(x+h)=f(x)+f’(x)*h+f’’(x)/2*h^2+f’’’(c1)/6*h^3

f(x-h)=f(x)-f’(x)*h+f’’(x)/2*h^2-f’’’(c1)/6*h^3

So: f(x+h)-f(x-h)=2*f’(x)*h+h^3/6*(f’’’(c1)+f’’’(c2))

Then, f’(x)(f(x+h)-f(x-h))/2h

(a)  f’(0)=5*(e^0.1-e^(-0.1))

(b) f’’(0)=50*(e^0.01-e^(-0.01))

(c)  f’’’(0)=500*(e^0.001-e^(-0.001))

Computer Problems

1.     Make a table of the error of the three-point centered-difference formula for f’(0), where f(x) = sinx − cosx, with h = 10 −1 ,...,10 −12 , as in the table in Section 5.1.2. Draw a plot of the results. Does the minimum error correspond to the theoretical expectation?

Solution:

Because the three-point centered-difference formula is            f’(x)(f(x+h)-f(x-h))/2h=[sin(x+h)-sin(x-h)+cos(x-h)-cos(x+h)]/2h

f’(0) (sin(h)-sin(-h)+cos(-h)-cos(h))/2*h

=(2*sin(h))/2*h

=sin(h)/h

f’(0)=cos(0)+sin(0)=1

Code:

m=-1:-1:-12

h(-m)=10.^m

f(-m)=sin(h)/h

 h formula error 10^(-1) 0.998334166468282 0.001665833531718 10^(-2) 0.999983333416667 0.000016666583333 10^(-3) 0.999999833333342 0.000000166666658 10^(-4) 0.999999998333333 0. 000000001666667 10^(-5) 0.999999999983333 0. 000000000016667 10^(-6) 0.999999999999833 0. 000000000000167 10^(-7) 0.999999999999998 0. 000000000000002 10^(-8) 1 0 10^(-9) 1 0 10^(-10) 1 0 10^(-11) 1 0 10^(-12) 1 0

Chap 5.2

Exercises

1.     Apply the composite Trapezoid Rule with m = 1,2, and 4 panels to approximate the integral.Compute the error by comparing with the exact value from calculus.

Solution:

(b) When m=1, = Π/4-Π^3/8/12*(-cos(c))

[0.462416114644325, 0.785398163397448]

When m=2, = Π/8*(1+0)-Π^3/32/12*(-cos(c))[0.311953569510443,0.392699081698724]

When m=4, = Π/16*(1+2*(cos(pi/8)+ cos(2*pi/8)+ cos(3*pi/8))+0)-Π^3/128/12*(-cos(c))[0.966929422925705, 0.987115800972776]

But,-sin(pi/2)+sin(0)=-1

2.      Apply the Composite Midpoint Rule with m = 1,2, and 4 panels to approximate the integrals in Exercise 1, and report the errors.

Solution:

(b)When m=1,=pi/2*(1/2)^(1/2)+pi^3/8/24*cos(c)[1.110720734539592, 1.272211758916153]

When m=2,=pi/4*(cos(pi/6)+cos(pi/3))+pi^3/32/24*cos(c)[1.072873843286556, 1.113246599380696]

When m=4,=pi/8*(cos(pi/10)+cos(pi/5)+ cos(3*pi/10)+ cos(2*pi/5))+pi^3/128/24*cos(c)[1.043352670094047, 1.053445859117582]

3.     Apply the composite Simpson’s Rule with m = 1,2, and 4 panels to the integrals in Exercise 1,and report the errors.

Solution:

(b) When m=1,h=pi/4,=1/3*(cos(0)+cos(pi/2)+4*cos(pi/4))+pi^5/256/360*sin(c)(1.276142374915397,1.279462901008987)

When

m=2,h=pi/8,=1/3*(cos(0)+cos(pi/2)+4*(cos(pi/8)+cos(3*pi/8))+2*cos(pi/4))+pi^5/64/64/360*sin(c)(2.546821807292867, 2.547029340173716)

When

m=4,h=pi/16,=1/3*(cos(0)+cos(pi/2)+4*(cos(3*pi/16)+cos(5*pi/16)+cos(7*pi/16))+2*(cos(pi/8)+cos(pi/4)+cos(6*pi/16)))+pi^5/256/256/360*sin(c)(3.785286720492984, 3.785299691298037)

Computer Problems

1.     Use the composite Trapezoid Rule with m = 16 and 32 panels to approximate the definite integral. Compare with the correct integral and report the two errors.

Solution:

Code:

function int_f = ComTrapezium(f,a,b,m)

% input :   f : function handler

%           a : the lower limit of integral

%           b : the upper limit of integral

%           m : cut integral area into m peace

% output :  int_f : the answer of the integral

h=(b-a)/m;

int_f=0;

if m>=2

for i=1:m-1

int_f=int_f+2*f(a+h*i);

end

end

int_f=int_f+f(a)+f(b);

int_f=int_f*h/2;

end

(d) the correct integral is:

>> syms x;

>> f=(x.^2).*(log(x));

>> int(f,1,3)

ans =

9*log(3) - 26/9

M=16

>> f = @(x) (x.^2).*(log(x));

>> int4_16 = ComTrapezium(f, 1, 3, 16)

int4_16 =

7.009809235924680

Error(16)= 9*log(3) - 26/9-7.009809235924680=-0.011187526800584

M=32

>> f = @(x) (x.^2).*(log(x));

>> int4_32 = ComTrapezium(f, 1, 3, 32)

int4_32 =

7.001418506166504

Error(32)= 9*log(3) - 26/9-7.001418506166504=-0.002796797042408

(h) the correct integral is:

>> syms x;

>> f = x./((x.^4+1).^0.5);

>> int(f,0,1)

ans =

log(2^(1/2) + 1)/2

M=16

>> f = @(x) x./((x.^4+1).^0.5);

>> int8_16=ComTrapezium(f, 0, 1, 16)

int8_16 =

0.440361182629694

Error(16)= log(2^(1/2) + 1)/2-0.440361182629694= 3.256108800774871e-04

M=32

>> f = @(x) x./((x.^4+1).^0.5);

>> int8_32=ComTrapezium(f, 0, 1, 32)

int8_32=

0.440605407679783

Error(32)= log(2^(1/2) + 1)/2-0.440605407679783= 8.138582998845623e-05

2.     Apply the composite Simpson’s Rule to the integrals in Computer Problem 1. Use m = 16 and 32, and report errors.

Solution:

(d)Code:

function y=f(x)

y=(x.^2).*(log(x));

function S_n=S_P_S(a,b,n)

h=(b-a)/n;

for k=0:n

x(k+1)=a+k*h;

x_k(k+1)=x(k+1)+1/2*h;

if (x(k+1)==0)||(x_k(k+1)==0)

x(k+1)=10^(-10);

x_k(k+1)=10^(-10);

end

end

S_1=h/6*(f(x(1))+f(x(n+1)));

for i=2:n

F_1(i)=h/3*f(x(i));

end

for j=1:n

F_2(j)=2*h/3*f(x_k(j));

end

S_2=sum(F_1)+sum(F_2);

S_n=S_1+S_2;

M=16

>> S_n=S_P_S(1,3,16)

S_n =

6.998621596247110

Error(16)=9*log(3)-26/9-6.998621596247110=1.128769859803924e-07

M=32

>> S_n=S_P_S(1,3,32)

S_n =

6.998621702062213

Error(32)= 9*log(3)-26/9-6.998621702062213=7.061883522396784e-09

(h) Code:

function y=f(x)

y= x./((x.^4+1).^0.5);

function S_n=S_P_S(a,b,n)

h=(b-a)/n;

for k=0:n

x(k+1)=a+k*h;

x_k(k+1)=x(k+1)+1/2*h;

if (x(k+1)==0)||(x_k(k+1)==0)

x(k+1)=10^(-10);

x_k(k+1)=10^(-10);

end

end

M=16

>> S_n=S_P_S(0,1,16)

S_n =

0.439384733322569

Error(16)= log(2^(1/2) + 1)/2-0.439384733322569= 0.001302060187202

M=32

>> S_n=S_P_S(0,1,32)

S_n =

0.440361274094287

Error(32)=log(2^(1/2)+1)/2-0.440361274094287=3.255194154844765e-04

3.     Use the composite Trapezoid Rule with m = 16 and 32 panels to approximate the definite integral.

Solution:

Code:

function int_f = ComTrapezium(f,a,b,m)

% input :   f : function handler

%           a : the lower limit of integral

%           b : the upper limit of integral

%           m : cut integral area into m peace

% output :  int_f : the answer of the integral

h=(b-a)/m;

int_f=0;

if m>=2

for i=1:m-1

int_f=int_f+2*f(a+h*i);

end

end

int_f=int_f+f(a)+f(b);

int_f=int_f*h/2;

end

(a)m=16

>> f = @(x)(eps(x.^2));

>> int_f = ComTrapezium(f,0,1,16)

int_f =

5.448115916739660e-17

m=32

>> f = @(x)(eps(x.^2));

>> int_f = ComTrapezium(f,0,1,32)

int_f =

5.278031700930996e-17

(f)m=16

>> f = @(x)(cos(eps(x)));

>> int_f = ComTrapezium(f,0,pi,16)

int_f =

3.141592653589793

m=32

>> f = @(x)(cos(eps(x)));

>> int_f = ComTrapezium(f,0,pi,32)

int_f =

3.141592653589793

4.     Apply the composite Simpson’s Rule to the integrals of Computer Problem 3, using m = 16 and 32.

Solution:

(a)Code:

function y=f(x)

y= eps(x.^2);

function S_n=S_P_S(a,b,n)

h=(b-a)/n;

for k=0:n

x(k+1)=a+k*h;

x_k(k+1)=x(k+1)+1/2*h;

if (x(k+1)==0)||(x_k(k+1)==0)

x(k+1)=10^(-10);

x_k(k+1)=10^(-10);

end

end

m=16

>> S_n=S_P_S(0,1,16)

S_n =

5.220433460517703e-17

m=32

>> S_n=S_P_S(0,1,32)

S_n =

5.076099046305572e-17

(h)Code:

function y=f(x)

y=cos(eps(x));

function S_n=S_P_S(a,b,n)

h=(b-a)/n;

for k=0:n

x(k+1)=a+k*h;

x_k(k+1)=x(k+1)+1/2*h;

if (x(k+1)==0)||(x_k(k+1)==0)

x(k+1)=10^(-10);

x_k(k+1)=10^(-10);

end

end

S_1=h/6*(f(x(1))+f(x(n+1)));

for i=2:n

F_1(i)=h/3*f(x(i));

end

for j=1:n

F_2(j)=2*h/3*f(x_k(j));

end

S_2=sum(F_1)+sum(F_2);

S_n=S_1+S_2;

m=16

>> S_n=S_P_S(0,pi,16)

S_n =

3.141592653589794

m=32

>> S_n=S_P_S(0,pi,32)

S_n =

3.141592653589792