﻿ 17074217数值分析9 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 17074217数值分析9

5.1 Exercises

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

(a)

The three-point centered-difference formula evaluates to

f(x)=[f(x+h)-f(x-h)]/2h=[e^0.1-e^(-0.1)]/0.2= 1.001667500198441.

(b)

The three-point centered-difference formula evaluates to

f(x)=[f(x+h)-f(x-h)]/2h=[e^0.01-e^(-0.01)]/0.02= 1.000016666750000.

(c)

The three-point centered-difference formula evaluates to

f(x)=[f(x+h)-f(x-h)]/2h=[e^0.001-e^(-0.001)]/0.002=  1.000000166666681.

5.1 Computer Problems

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

Code:

function E=tcf(f,x)

for i=1:12

h(i)=1/(10.^i);

end

for i=1:12

F(i)=(f(x+h(i))-f(x-h(i)))/(2*h(i));

e(i)=F(i)-1;

end

E(1:12,1)=h(1:12);

E(1:12,2)=F(1:12);

E(1:12,3)=e(1:12);

E=vpa(E,14);

x0=h(1:12);

y0=e(1:12);

plot(x0,y0)

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

E=tcf(f,0)

E =

[            0.1, 0.99833416646828,         -0.0016658335317179]

[           0.01, 0.99998333341667,       -0.000016666583334768]

[          0.001, 0.99999983333338,      -0.0000001666666240574]

[         0.0001, 0.99999999833289,   -0.0000000016671100055987]

[        0.00001, 0.99999999998435, -0.000000000015653367491097]

[       0.000001, 0.99999999997324, -0.000000000026755486715047]

[      0.0000001, 0.99999999947364,  -0.00000000052635584779637]

[     0.00000001, 0.99999999947364,  -0.00000000052635584779637]

[    0.000000001,  1.0000000272292,     0.000000027229219767833]

[   0.0000000001,  1.0000000827404,      0.00000008274037099909]

[  0.00000000001,  1.0000000827404,      0.00000008274037099909]

[ 0.000000000001,  1.0000333894311,        0.000033389431109754]

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.

π/2
(b)
∫  cos(x)dx
0

m=1:
h=(π/2-0)/1=π/2 ,

π/2
∫ cos(x)dx=h/2*(y0+y1)=π/4*(cos(0)+cos(π/2))=0.7854,
0
error=(π/2-0)/12*(π/2).^2*|-cos(c)|≤ 1/12*(π/2).^3=0.3230,

π/2
∫ cos(x)dx=0.7854±0.3230;
0

m=2:
h=(π/2-0)/2=π/4 ,

π/2
∫ cos(x)dx=h/2*(y0+y2+2*y1)=π/8*(cos(0)+cos(π/2)+2*cos(π/4))=0.9481,
0
error=(π/2-0)/12*(π/4).^2*|-cos(c)|≤ (π/24)*(π.^2/16)=0.0807,

π/2
∫ cos(x)dx=0.9481±0.0807;
0

m=4:
h=(π/2-0)/4=π/8,

π/2
∫ cos(x)dx=h/2*(y0+y4+2*y1+2*y2+2*y3)=π/16*(cos(0)+cos(π/2)+2*cos(π/8)+2*cos(π/4)+2*cos(3*π/8))=0.9871,
0
error=(π/2-0)/12*(π/8).^2*|-cos(c)|≤ (π/24)*(π.^2/64)=0.0202,

π/2
∫ cos(x)dx=0.9871±0.0202.
0

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

in Exercise 1, and report the errors.

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.

m=1: h=(π/2-0)/(2*1)=π/4,

π/2
∫ cos(x)dx=h/3*(y0+y2+4*y1)=π/12*(cos(0)+cos(π/2)+4*cos(π/4))=1.0023,
0
error=(π/2-0)/180*(π/4).^4*|cos(c)|≤(π.^5)/(360*4.^4)=0.0033,

π/2
∫ cos(x)dx=1.0023±0.0033;
0

m=2: h=(π/2-0)/(2*2)=π/8,

π/2
∫ cos(x)dx=h/3*(y0+y4+4*y1+4*y3+2*y2)=π/24*(cos(0)+cos(π/2)+4*cos(π/8)+4*cos(3*π/8)+2*cos(π/4))=1.00013458,
0
error=(π/2-0)/180*(π/8).^4*|cos(c)|≤(π.^5)/(360*8.^4)=0.0002075329,

π/2
∫ cos(x)dx=1.00013458±0.0002075329;
0

m=4: h=(π/2-0)/(2*4)=π/16,

π/2
∫ cos(x)dx=h/3*(y0+y8+4*y1+4*y3+4*y5+4*y7+2*y2+2*y4+2*y6)
0         =π/48*(cos(0)+cos(π/2)+4*cos(π/16)+4*cos(3*π/16)+4*cos(5*π/16)+4*cos(7*π/16)+2*cos(π/8)+2*cos(π/4)+2*cos(3*π/8))
=1.0000083,
error=(π/2-0)/180*(π/16).^4*|cos(c)|≤(π.^5)/(360*16.^4)=1.297081*10.^(-5),

π/2
∫ cos(x)dx=1.0000083±1.297081*10.^(-5).
0

5.2 Computer Problems

1. Use the composite Trapezoid Rule with m = 16 and 32 panels to approximate the defifinite

integral. Compare with the correct integral and report the two errors.

3           1
(d):∫ x.^2*ln(x)dx;    (h):∫ x/(√(x.^4+1))dx.
1           0

Code

function y=TR(f,a,b,m)
h=(b-a)/m;
for i=1:m+1
x(i)=a+(i-1)*h;
end
y=f(x(1))+f(x(m+1));
for i=2:m
y=y+2*f(x(i));
end
y=y*h/2;

(d):

>>f=@(x) x.^2*log(x);
>> y=TR(f,1,3,16)

y =

7.0098

>>f=@(x) x.^2*log(x);
>> y=TR(f,1,3,32)

y =

7.0014

m=16:
error=(2/12)*(2/16).^2*|2*ln(c)+3|≤(1/(24*16))*(2*ln(3)+3)=0.0135;
m=32:
error=(2/12)*(2/32).^2*|2*ln(c)+3|≤(2/12)*(2/32).^2*(2*ln(3)+3)=0.0034;

(h):

>>f=@(x) x/((x.^4+1).^(1/2));
>> y=TR(f,0,1,16)

y =

0.4404

>>f=@(x) x/((x.^4+1).^(1/2));
>> y=TR(f,0,1,32)

y =

0.4406

m=16:
error=(1/12)*(1/16).^2*|(2*c.^7-10*c.^3)/(c.^4+1).^(5/2)|≤(1/(12*16))*(8/(2.^(5/2)))=0.00046;
m=32:
error=(1/12)*(1/32).^2*|(2*c.^7-10*c.^3)/(c.^4+1).^(5/2)|≤(1/12)*(1/32).^2*(8/(2.^(5/2)))=0.00012;

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

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.
1          π
(a):∫ e.^(x.^2)dx;     (f):∫ cos(e.^(x))dx.
0           0

Code:

function y=TR(f,a,b,m)
h=(b-a)/m;
for i=1:m+1
x(i)=a+(i-1)*h;
end
y=f(x(1))+f(x(m+1));
for i=2:m
y=y+2*f(x(i));
end
y=y*h/2;

Results:

(a):

>>f=@(x) exp(x.^2);
>> y=TR(f,0,1,16)

y =

1.4644

>>f=@(x) exp(x.^2);
>> y=TR(f,0,1,32)

y =

1.4631

(f):

>>f=@(x) cos(exp(x));
>> y=TR(f,0,pi,16)

y =

-0.2780

>>f=@(x) cos(exp(x));
>> y=TR(f,0,pi,32)

y =

-0.3568

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

and 32.

(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

6.1 Exercises

3. Use separation of variables to fifind solutions of the IVP given by y(0) = 1 and the following

differential equations:

(a) y′ = t (b) y′ = t2y (c) y′ = 2(t + 1)y(d) y′ = 5t4y (e) y′ = 1/y2(f) y′ = t3/y

(a)  dy/dt=t;

;

then  y=t2/2+c;

because y(0)=1,then c=1;

so y= t2/2+1

(b)  dy/dt=t2y;

;

then  ln|y|=t2/2+c;

because y(0)=1,then c=0;

so y=e t^2/2

(c)  dy/dt=2(t+1)y;

;

then  ln|y|=t2+2t+c;

because y(0)=1,then c=0;

so y=e t^2+2t

(d)  dy/dt=5t 4 y;

;

then  ln|y|=t5+c;

because y(0)=1,then c=0;

so y=e t^5

(e)  dy/dt=1/y 2;

;

then  y3/3=t+c;

because y(0)=1,then c=1/3;

so t=(y3-1)/3

(f)    dy/dt= t 3 /y 2;

;

then  y3/3=t4/4+c;

because y(0)=1,then c=1/3;

so t4/4= (y3-1)/3

5. Apply Euler’s Method with step size h = 1/4 to the IVPs in Exercise 3 on the interval

[0,1]. List the wi, i = 0,... ,4, and fifind the error at t = 1 by comparing with the correct

Solution

(a)  y ‘= t, the correct solution is y= t2/2+1,so

 wi W0 W1 W2 W3 W4 The correct wi 1 1.5 3 5.5 9

Use the euler method,

When i=0, w1=w0+1/4*0=1

When i=1, w2=w1+1/4*1=1.25

When i=2, w3=w2+1/4*2=1.75

When i=3, w4=w3+1/4*3=2.5

 wi W0 W1 W2 W3 W4 The euler wi 1 1 1.25 1.75 2.5

When i=1, the error is 0.5.

(b)  y ‘= t2y, the correct solution is y= e t^2/2,so

 wi W0 W1 W2 W3 W4 The correct wi 1 e1/2 e2 e9/2 e8

Use the euler method,

When i=0, w1=w0+1/4*0=1

When i=1, w2=w1+1/4*1*1=1.25

When i=2, w3=w2+1/4*4*1.25=2.5

When i=3, w4=w3+1/4*9*2.5= 8.125

 wi W0 W1 W2 W3 W4 The euler wi 1 1 1.25 2.5 8.125

When i=1, the error is 0.648721270700128.

(c)  y ‘= 2(t+1)y, the correct solution is y= e t^2+2t,so

 wi W0 W1 W2 W3 W4 The correct wi 1 e3 e8 e15 e24

Use the euler method,

When i=0, w1=w0+1/4*2=1.5

When i=1, w2=w1+1/4*4*1.5=3

When i=2, w3=w2+1/4*6*3=7.5

When i=3, w4=w3+1/4*8*7.5=22.5

 wi W0 W1 W2 W3 W4 The euler wi 1 1.5 3 7.5 22.5

When i=1, the error is 17.085536923187668

(d)  y ‘= 5t 4 y, the correct solution is y= e t^5,so

 wi W0 W1 W2 W3 W4 The correct wi 1 e e32 e3^5 e4^5

Use the euler method,

When i=0, w1=w0+1/4*0=1

When i=1, w2=w1+1/4*5*1=2.25

When i=2, w3=w2+1/4*5*2^4*2.25=47.25

When i=3, w4=w3+1/4*5*3^4*47.25= 4.831312500000000e+03

 wi W0 W1 W2 W3 W4 The euler wi 1 1 2.25 47.25 4.831312500000000e+03

When i=1, the error is 1.718281828459046

(e)  y ‘= 1/y 2, the correct solution is t=(y3-1)/3,so

 wi W0 W1 W2 W3 W4 The correct wi 1 41/3 71/3 101/3 131/3

Use the euler method,

When i=0, w1=w0+1/4*1=1.25

When i=1, w2=w1+1/4*1/(1.25^2)= 1.410000000000000

When i=2, w3=w2+1/4*1/(1.41^2)= 1.535748201800714

When i=3, w4=w3+1/4*1/( 1.535748201800714^2)= 1.641746764812362

 wi W0 W1 W2 W3 W4 The euler wi 1 1.25 1.41 1.535748201800714 1.641746764812362

When i=1, the error is 0.337401051968199

(f)    y ‘= t 3 /y 2, the correct solution is t4/4= (y3-1)/3,so

 wi W0 W1 W2 W3 W4 The correct wi 1 (7/4)1/3 131/3 (247/4)1/3 491/3

Use the euler method,

When i=0, w1=w0+1/4*1=1.25

When i=1, w2=w1+1/4*1/(1.25^2)= 1.410000000000000

When i=2, w3=w2+1/4*8/(1.41^2)= 2.415985614405714

When i=3, w4=w3+1/4*27/(2.415985614405714^2)= 3.572404258290336

 wi W0 W1 W2 W3 W4 The euler wi 1 1.25 1.41 2.415985614405714 3.572404258290336

When i=1, the error is -0.044928867912385

6.2 Exercises

1. Using initial condition y(0) = 1 and step size h = 1/4, calculate the Trapezoid Method

approximation w0,... ,w4 on the interval [0,1]. Find the error at t = 1 by comparing with the

correct solution found in Exercise 6.1.3

(a) y′ = t (b) y′ = t2y (c) y′ = 2(t + 1)y

(a)the correct solution is y=t2/2+1

 wi W0 W1 W2 W3 W4 Correct answer 1 1.5 3 5.5 9

Use the Trapezoid Method wi+1=wi+h/2(f(ti,yi)+f(ti+h,wi+hf(ti,yi)))

When i=0, w1=w0+1/8*(0+1/4)= 1.03125

When i=1, w2=w1+1/8*(1/4+1/2)= 1.125

When i=2, w3=w2+1/8*(1/2+1/2+1/4)= 1.28125

When i=3, w4=w3+1/8*(3/4+3/4+1/4)=1.5

 wi W0 W1 W2 W3 W4 Trapezoid Method answer 1 1.03125 1.125 1.28125 1.5

When i=1,the error is 0.46875

(a)  the correct solution is y=et^3/3

 wi W0 W1 W2 W3 W4 Correct answer 1 e1/3 e8/3 e9 e64/3

Use the Trapezoid Method wi+1=wi+h/2(f(ti,yi)+f(ti+h,wi+hf(ti,yi)))

When i=0, w1=w0+1/8*(0+1/16)= 1. 0078125

When i=1, w2=w1+1/8*(1/16*1.0078125+1/4*1.0078125*65/64)= 1.047672271728516

When i=2, w3=w2+1/8*(1/4*1.047672271728516+9*17/256*1.047672271728516)= 1.158680515363813

When i=3, w4=w3+1/8*(9/16*1.158680515363813+73/64*1.158680515363813)= 1.405352734454937

 wi W0 W1 W2 W3 W4 Trapezoid Method answer 1 1.0078125 1.047672271728516 1.158680515363813 1.405352734454937

When i=1,the error is 0.387799925086090

(b)  the correct solution is y=et^2+2t

 wi W0 W1 W2 W3 W4 Correct answer 1 e3 e8 e15 e24

Use the Trapezoid Method wi+1=wi+h/2(f(ti,yi)+f(ti+h,wi+hf(ti,yi)))

When i=0, w1=w0+1/8*(2+15/4)= 1. 71875

When i=1, w2=w1+1/8*(4.21875+3*(1.71875+1/4*4.21875))= 3.2861328125

When i=2, w3=w2+1/8*(3*3.2861328125+49/8*3.2861328125)= 7.034378051757813

When i=3, w4=w3+1/8*(11*7.034378051757813)= 16.706647872924810

 wi W0 W1 W2 W3 W4 Trapezoid Method answer 1 1. 71875 3.2861328125 7.034378051757813 16.706647872924810454937

When i=1,the error is 18.366786923187668