﻿ 17074188 赵文轩 作业9 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 17074188 赵文轩 作业9

P252

5.1Exercises

1.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.

(a)h=0.1,f'(0)(e^0.1-e^(-0.1))/0.2=5*(e^0.1-e^(-0.1)),f'(x)=f''(x)=f'''(x)=e^x,-0.1<c<0.1,f'''(c)=e^0=1,error=(-h^2*f'''(c))/6 =(-0.1^2*'''(c))/6=-0.0017.

(b)h=0.01,f'(0)(e^0.01-e^(-0.01))/0.02=50*(e^0.01-e^(-0.01)),error=(-h^2*f'''(c))/6 =(-0.01^2*f'''(c))/6=-1.67*10^(-5).

(c)h=0.001,f'(0)(e^0.001-e^(-0.001))/0.002=500*(e^0.001-e^(-0.001)),error=(-h^2*f'''(c))/6=(-0.001^2*f'''(c))/6=-1.67*10^(-7).

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?

(1)Code:

function TPCDF

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

A=zeros(12,3);

for i=1:1:12

h=10^(-1*i);

y1=(f(0+h)-f(0-h))/(2*h);

y2=1-y1;

A(i,1)=h;

A(i,2)=y1;

A(i,3)=y2;

end

disp(A);

>> x=[10^(-1);10^(-2);10^(-3);10^(-4);10^(-5);10^(-6);10^(-7);10^(-8);10^(-9);10^(-10);10^(-11);10^(-12) ];

y=[0.001665833531718;0.000016666583335; 0.000000166666624; 0.000000001667110;0.000000000015653; 0.000000000026755;0.000000000526356;0.000000000526356;-0.000000027229220;-0.000000082740371;-0.000000082740371;-0.000033389431110];

>> plot(x,y)

Result>> TPCDF

0.100000000000000   0.998334166468282   0.001665833531718

0.010000000000000   0.999983333416665   0.000016666583335

0.001000000000000   0.999999833333376   0.000000166666624

0.000100000000000   0.999999998332890   0.000000001667110

0.000010000000000   0.999999999984347   0.000000000015653

0.000001000000000   0.999999999973245   0.000000000026755

0.000000100000000   0.999999999473644   0.000000000526356

0.000000010000000   0.999999999473644   0.000000000526356

0.000000001000000   1.000000027229220  -0.000000027229220

0.000000000100000   1.000000082740371  -0.000000082740371

0.000000000010000   1.000000082740371  -0.000000082740371

0.000000000001000   1.000033389431110  -0.000033389431110

P263

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.

Y=cos(x);

y0=cos0=1, y1=cos1, y2=cos2, y4=cos4;

(a) m=1

f0 π/2 cosxdx=(h/2)*(y0+y1+2y0)-π/24*h^2f”(c)=0.785398;

Error=0.214602;

(a) m=2

f0 π/2 cosxdx=(h/2)*(y0+y2+2(y0+y1))-π/24*h^2f”(c)=0.948059;

error=0.051941;

(a) m=4

f0 π/2 cosxdx=(h/2)*(y0+y4+2(y0+y1+y2+y3))-π/24*h^2f”(c)=0.987116;

error=0.012884.

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

fx0 x2 f(x)dx=h/3(y0+4y1+y2)-h^5/90*f^(4)(c)

(a) m=1

fx0 x2 f(x)dx=h/3(y0+4y1+y2)-h^5/90*f^(4)(c)=1.002280

Error=0.002280;

(a) m=2

f0 π/2 cosxdx=(h/2)*(y0+y2+2(y0+y1))-π/24*h^2f”(c)=1.000135;

error=0.000135;

(a) m=4

f0 π/2 cosxdx=(h/2)*(y0+y4+2(y0+y1+y2+y3))-π/24*h^2f”(c)=1.000008;

error=0.000008.

5.2computer 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.

(d) f1 3 x^2*lnx dx

(h) f0 1 x/(x^4+1)^1/2 dx

(d) Code

function t=CT(a,b,h)

format long

x=a:h:b;

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

y(1)=0;

t=0;

for k=1:(b-a)/h,

t= t+y(k)+y(k+1);

end

t=t*h/2;

>> CT(1,3,0.1)

ans =

7.005781622292308

(h)

function t=CT(a,b,h)

format long

x=a:h:b;

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

y(1)=0;

t=0;

for k=1:(b-a)/h,

t= t+y(k)+y(k+1);

end

t=t*h/2;

>>  CT(0,1,0.1)

ans =

0.195723914362966

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

(a)f0 1 e^x^2 dx

(f)f0 π cose^x dx

(a)

function t=CT(a,b,h)

format long

x=a:h:b;

y=exp(x.^2) ;

y(1)=0;

t=0;

for k=1:(b-a)/h,

t= t+y(k)+y(k+1);

end

t=t*h/2;

>> CT(0,1,0.1)

ans =

1.417174692738799

(f)

function t=CT(a,b,h)

format long

x=a:h:b;

y=cos(exp(x)) ;

y(1)=0;

t=0;

for k=1:(b-a)/h,

t= t+y(k)+y(k+1);

end

t=t*h/2;

>>  CT(0,pi,0.1)

ans =

-0.367124048484025

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

differential equations:

(a) y  = t                                     (d) y  = 5t 4 y

Solution：（aAssuming that y=!0,divide both sides by y,sparate variables,and integrate ,

as follows

y=ct+k;

(b) Assuming that y=!0,divide both sides by y,sparate variables,and integrate ,

as follows

dy/y=5t

ln|y|=5/2t^2+k

y=exp(5/2t^2+k)=exp(5/2t^2)exp(k)

the initial condition y(0)=1 implies y=exp(5/2t^2).