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

17074211作业9

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.

f’(x)=(f(x+h)-f(x-h))/2h=(e^(x+h)-e^(x-h))/2h;

(a)f’(0)=(e^0.1-e^(-0.1))/2*0.11.007

(b)f’(0)=(e^0.01-e^(-0.01))/2*0.011.0000167

(c) f’(0)=(e^0.001-e^(-0.001))/2*0.0011.0000002.

5.1 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 5.1.2.Draw a plot of the results.Does the minimum error correspond to the theoretical expectation?

(1)Code:

>> x=0;

for n=-12:1:-1

h=10^(n)

f=vpa((sin(x+h)-cos(x+h)-sin(x-h)+cos(x-h))/(2*h),11)

error=vpa(1-f,11)

end

h =

1.0000e-012

f =

.99997787828

error =

.2212172e-4

h =

1.0000e-011

f =

1.0000000827

error =

-.827e-7

h =

1.0000e-010

f =

1.0000000827

error =

-.827e-7

h =

1.0000e-009

f =

.99999997172

error =

.2828e-7

h =

1.0000e-008

f =

1.0000000050

error =

-.50e-8

h =

1.0000e-007

f =

.99999999947

error =

.53e-9

h =

1.0000e-006

f =

1.0000000000

error =

0.

h =

1.0000e-005

f =

.99999999998

error =

.2e-10

h =

1.0000e-004

f =

.99999999833

error =

.167e-8

h =

1.0000e-003

f =

.99999983333

error =

.16667e-6

h =

0.0100

f =

.99998333342

error =

.1666658e-4

h =

0.1000

f =

.99833416647

error =

.166583353e-2

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

>>y0=[.2212172e-4 -.827e-7 -.827e-7 .2828e-7 -.50e-8 .53e-9 0 .2e-10 .167e-8 .16667e-6 .1666658e-4 .166583353e-2];

h                              f                                error

10^(-1)              .99833416647            .166583353e-2

10^(-2)              .99998333342            .1666658e-4

10^(-3)              .99999983333            .16667e-6

10^(-4)              .99999999833            .167e-8

10^(-5)              .99999999998            .2e-10

10^(-6)              1.0000000000            0.

10^(-7)              .99999999947            .53e-9

10^(-8)              1.0000000050            -.50e-8

10^(-9)              .99999997172            .2828e-7

10^(-10)            1.0000000827             -.827e-7

10^(-11)            1.0000000827             -.827e-7

10^(-12)            .99997787828             .2212172e-4

>> c=newtdd(x0,y0,12);

>> x=10^(-12):10^(-1);

>>  y=nest(11,c,x,x0);

>> plot(x0,y0,'o',x,y)

(2)Result:

error=-(h^2/6)(f’’’(c))

The minimum error closes to the theoretical expectation.

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.

(b)(cosx)(0,π/2)

(cosx)(0,π/2)=1-1=0.

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

(cosx)(0,π/2)=(π/2)/2(y0+y1)=π/4(1+0)=π/4

error=1-π/40.2146

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

(cosx)(0,π/2)=(π/4)/2(y0+y2+2y1)=π/8(1+0+2^(1/2))0.9481

error=1-0.9481=0.0519

(3)m=4, h=(b-a)/m=(π/2-0)/4=π/8

(cosx)(0,π/2)=(π/8)/2(y0+y4+2(y1+y2+y3))

=π/16[1+0+2(cos(π/8)+ cos2π/8)+ cos(3π/8))]

0.9781

error=1-0.9781=0.0129.

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=1,h=(π/2-0)/1=π/2:

π/2)=sqrt(2)*pi/4;

So the error is 1- sqrt(2)*pi/4=-0.110720734539592;

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

π/4)= 1.0261;

So the error is 1- 1.0261=-0.0261;

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

π/8)= 1.0064;

So the error is 1- 1.0064=-0.0064;

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

(cosx)(0,π/2)=1-0=1

(1)m=1,h=(b-a)/2m=(π/2-0)/2=π/4

(cosx)(0,π/2)=(π/4)/3(y0+y2+4y1)=π/12(1+0+4*(2)^(1/2)/2)1.0022799

error=|1-1.0022799|=0.0022799

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

(cosx)(0,π/2)=(π/8)/3(y0+y4+4(y1+y3)+2y2)≈1.0001346

error=|1-1.0001346|=0.0001346.

(3)m=4, h=(b-a)/2m=(π/2-0)/8=π/16

(cosx)(0,π/2)=(π/16)/3(y0+y8+4(y1+y3+y5+y7)+2(y2+y4+y6))1.0000083

error=|1-1.0000083|=0.0000083.

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

(b):

When m=1,Simpson’s Rule sets h=(1-0)/2*1=1/2.The approximation is

1/2/3(y0+y2+4+2)= 0.7833

The exact value from calculus is  =arctanx|(0-1)=pi/4;

So the error is pi/4-0.7833= 0.002098163397448

When m=2,Simpson’s Rule sets h=(1-0)/2*2=1/4.The approximation is

1/4/3(y0+y4+4+2)= 0.785392156862745

So the error is pi/4-0.785392156862745= 6.006534703284494e-06

When m=4,Simpson’s Rule sets h=(1-0)/2*4=1/8.The approximation is

1/8/3(y0+y8+4+2)= 0.785398125614677

So the error is pi/4-0.785398125614677= 3.778277124499851e-08

5.2computer programs

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

(d):  (h)

Code:

function y=CTR(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);

>> a=1;b=3;m=16;

>> y=CTR(f,a,b,m)

Result

y =

7.0098

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

a=1;b=3;m=32;

y=CTR(f,a,b,m)

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));

>> a=0;b=1;m=16;

>> y=CTR(f,a,b,m)

Result

y =

0.4404

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

a=0;b=1;m=32;

y=CTR(f,a,b,m)

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 and 32, and report errors.

Code:

function y=CSR(f,a,b,m)

h=(b-a)/(2*m);

for i=1:2*m+1

x(i)=a+(i-1)*h;

end

y=f(x(1))+f(x(2*m+1));

for i=2:m

y=y+2*f(x(2*i-1));

end

for i=2:m+1

y=y+4*f(x(2*i-2));

end

y=y*(h/3);

Result:

(d):

m=16:

>> format long

>> y=CSR(f,1,3,16)

y =

6.998621596247110

m=32:

>> y=CSR(f,1,3,32)

y =

6.998621702062211

(h):

m=16:

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

>> y=CSR(f,0,1,16)

y =

0.440686816029813

m=32:

>> y=CSR(f,0,1,32)

y =

0.440686794915315

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

(a) (e^(x^2))(0,1)   (f)cos(e^x)(0,π)

The definite integral is 1.4644 and 1.4631.

(2)Code:

function trapezoid = trapezoid(a,b,m,f)

trapezoid=0;

for i=1:m-1

trapezoid=trapezoid+f(a+i*(b-a)/m);

end

trapezoid=(b-a)/(2*m)*(f(a)+f(b)+2*trapezoid);

end

(3)Result:

>> f=@(x)  exp(x^2);

>> trapezoid(0,1,16,f)

ans =

1.4644

>> trapezoid(0,1,32,f)

ans =

1.4631

The definite integral is -0.2780 and -0.3568.

(2)Code:

function trapezoid = trapezoid(a,b,m,f)

trapezoid=0;

for i=1:m-1

trapezoid=trapezoid+f(a+i*(b-a)/m);

end

trapezoid=(b-a)/(2*m)*(f(a)+f(b)+2*trapezoid);

end

(3)Result:

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

>> trapezoid(0,pi,16,f)

ans =

-0.2780

>> trapezoid(0,pi,32,f)

ans=

-0.3568

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

Code:

function y=CSR(f,a,b,m)

h=(b-a)/(2*m);

for i=1:2*m+1

x(i)=a+(i-1)*h;

end

y=f(x(1))+f(x(2*m+1));

for i=2:m

y=y+2*f(x(2*i-1));

end

for i=2:m+1

y=y+4*f(x(2*i-2));

end

y=y*(h/3);

Result:

(a):

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

>> y=CSR(f,0,1,16)

y =

1.462652033425411

>> y=CSR(f,0,1,32)

y =

1.462651763901493

(f):

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

>> y=CSR(f,0,pi,16)

y =

-0.383048399454781

>> y=CSR(f,0,pi,32)

y =

-0.376327835672488

5.3Exercises

1. Apply Romberg Integration to ﬁnd R33 for the integrals.

(c):

h1=1;h2=1/2;h2=1/4;h3=1/8;

R11=h1*(f(0)+f(1))/2=1.859140914229523;

R21=h2*(f(0)+f(1)+2(f(1/2))/2=1.753931092464826;

R31=R21/2+h3*(f(h3)+f(3*h3))= 1.727221904557517;

R22=(4R21-R11)/3=1.718861151876593;

R32=(4R31-R21)/3=1.718318841921747;

R33=(16R32-R22)/15=1.718282687924758

R=

1.859140914229523                   0                   0

1.753931092464826   1.718861151876593                   0

1.727221904557517   1.718318841921747   1.718282687924758

6.1 Exercises

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

(a)y'=t;      (b) y=t^2*y;        (c) y=2(t + 1)y

(d)y'=5*t^4*y;     (e) y=1/y^2 ;        (f) y=t^3/y^2

(a)

y'=t,dy/dt=t,y=t^2/2+c;

y(0)=1,c=1;

So y=t^2/2+1.

(b)

y=t^2*y, dy/dt=t^2*y, lny=t^3/3+c,y=ce^(t^3)

y(0)=1,c=1;

So |y|= e^(t^3) or y=0.

(c)

y=2(t + 1)y, dy/dt=2(t + 1)y, lny=t^2+2t+c, y=ce^(t^2+2t)

y(0)=1,c=1;

So  |y|=ce^(t^2+2t) or y=0.

(d)

y'=5*t^4*y,dy/dt=5*t^4*y(y!=0),ln(|y|)=t^5+c;

y(0)=1,c=0;

|y|=e^(t^5);

y=0,y'=0;

So |y|=e^(t^5) or y=0.

(e)

y=1/y^2,dy/dt=1/y^2, -1/y=t+c,

y(0)=1,c=-1;

So y=-1/(t-1)

(f)

y=t^3/y^2, dy/dt= t^3/y^2, -1/y=t^4/4+c

y(0)=1,c=-1;

So y=-1/(t^4/4-1)

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 ﬁnd the error at t =1 by comparing with the correct solution.

Euler’s Method

w0=y0=1

wi+1=wi+hf(ti,wi)

(a) w=[1.0000,1.0000,1.0625,1.1875,1.3750],error=0.1250

(b) w=[1.0000,1.0000,1.0156,1.0791,1.2309],error=0.1648

(c) w=[1.0000,1.5000,2.4375,4.2656,7.9980]error=12.0875

(d) w=[1.0000,1.0000,1.0049,1.0834,1.5119],error=1.2064

(e) w=[1.0000,1.2500,1.4100,1.5357,1.6417],error=0.0543

(f) w=[1.0000,1.0000,1.0039,1.0349,1.1334],error=0.0717

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=t^2y;        (c) y=2(t + 1)y

(d)y'=5*t^4*y;     (e) y=1/y^2 ;        (f) y=t^3/y^2

w0=y0=1

wi+1=wi+h/2(f(ti,wi)+f(ti+h,hf(ti,wi)))

(a) w=[1.0000,1.0313,1.1250,1.2813,1.5000],error=0

(b) w=[1.0000,1.0078,1.0477,1.1587,1.4054],error=0.0097

(c) w=[1.0000,1.7188,3.3032,7.0710,16.7935],error=3.2920

(d) w=[1.0000,1.0024,1.0442,1.3077,2.7068],error=0.0115

(e) w=[1.0000,1.2050,1.3570,1.4810,1.5817],error=0.0003

(f) w=[1.0000,1.0020,1.0193,1.0823,1.2182],error=0.0132