﻿ 17074192 作业九 Numerical Differentiation and Integration - 计算数学达人 - 专，学者，数值代数，微分方程数值解

## EXERCISES 5.1

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:

The three point centered-difference formula is f’(x) approximately equals (f(x+h)-f(x-h))/(2*h). then we can have

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

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

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

## COMPUTER PROBLEMS 5.1

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:

The code of the problem is

n=-linspace(1,12,12);

df0=ones(1,12);

h=10.^n;

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

for i=1:12

df(i)=(f(h(i))-f(-h(i)))/(2*h(i));

end

e=abs(df0-df);

plot(h,e);

The three point centered-difference formula is f’(x)(f(x+h)-f(x-h))/(2*h), and we can have that the exact f’(0)=1. The following table is given:

 h formula error 10^-1 0.998334166468282 0.00166583353171790 10^-2 0.999983333416665 1.66665833347679e-05 10^-3 0.999999833333376 1.66666624057399e-07 10^-4 0.999999998332890 1.66711000559872e-09 10^-5 0.999999999984347 1.56531454464925e-11 10^-6 0.999999999973245 2.67554867150466e-11 10^-7 0.999999999473644 5.26355847796367e-10 10^-8 0.999999999473644 5.26355847796367e-10 10^-9 1.00000002722922 2.72292197678325e-08 10^-10 1.00000008274037 8.27403709990904e-08 10^-11 1.00000008274037 8.27403709990904e-08 10^-12 1.00003338943111 3.33894311097538e-05

From the table, we can get that the minimum error doesn’t correspond to the theoretical expectation.

## EXERCISES 5.2

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.（cos(x)在[0,pi/2]上的积分）

Solution:

The rule of composite Trapezoid is

The exact value of the integral

is 1.

When m=1, we can have that h=pi/2, y0=1,ym=0; it means that pi*(0+1)/4=0.7854

When m=2, we can have that h=pi/4, y0=1,y1=sqrt(2)/2,y2=0; it means that pi*(1+0+sqrt(2))/8= 0.9481

When m=4,we can have that h=pi/8,y0=1,y1=cos(pi/8),y2=cos(pi/4),y3=cos(3*pi/8) and y4=0;it means that pi*(y0+y4+2*(y1+y2+y3))/16= 0.9871

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:

The composite Midpoint Rule is  where h=(b-a)/m and wi are the midpoints of the m equal subintervals of [a, b].

We have known that the exact value of the integral is 1.

When m=1, we can have that h=pi/2,f(w1)=sqrt(2)/2; it means that pi*sqrt(2)/4= 1.1107. And the error is 0.1107.

When m=2, we can have that h=pi/4, f(w1)=cos(pi/8), f(w2)=cos(3*pi/8); it means that pi*(cos(pi/8)+cos(3*pi/8))/4= 1.0262.  The error is 0.0262.

When m=4,we can have that h=pi/8, f(w1)=cos(pi/16),f(w2)=cos(3*pi/16),f(w3)=cos(5*pi/16) and f(w4)=cos(7*pi/16); it means that pi*( cos(pi/16)+cos(3*pi/16)+cos(5*pi/16)+ cos(7*pi/16))/8=1.0065.   The error is 0.0065.

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: the composite Simpson rule is

We have known that the exact value of the integral is 1.

When n=1, we can have that h=pi/4, y0=1, y1=sqrt(2)/2, y2=0; it means that pi*(1+0+2*sqrt(2))/12=1.0023. The error is 0.0023

When n=2, we can have that h=pi/8, y0=1,y1=cos(pi/8),y2=cos(pi/4),y3=cos(3*pi/8) and y4=0; it means that pi*(1+0+4*(cos(pi/8)+cos(3*pi/8))+2*(cos(pi/4)))/24=1.0001. The error is 0.0001.

When n=4, we can have that h=pi/16, y0=1, y1=cos(pi/16), y2=cos(pi/8), y3=cos(3*pi/16),y4=cos(pi/4),y5=cos(5*pi/16), y6=cos(3*pi/8),y7=cos(7*pi/16) and y8=1; then we can get that pi*(1+0+4*(y1+y3+y5+y7)+2*(y2+y4+y6))/48= 1.0000083. The error is 0. 0000083.

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

Solution: the composite Simpson rule is

And we have known that the exact value of the integral is arctan(1)=pi/40.785398163397448

When n=1, we can have that h=1/2, y0=1,y1=4/5 and y2=1/2; then we can get that 0.7833; the error is 0.0021.

When n=2, we can have that h=1/4, y0=1,y1=16/17,y2=4/5,y3=16/25 and y4=1/2;then we can have that (1+0.5+4*(y1+y3)+2*(y2))/12=0.785392156862745.error is 6.0065e-06

When n=4, we can have that h=1/8,y0=1,y1=64/65,y2=16/17,y3=64/73,y4=4/5,y5=64/89, y6=16/25, y7=64/113 and y8=1/2.; then we can have that (1+1/2+4*(y1+y3+y5+y7)+2*(y2+y4+y6))/24=0.785398125614677; the error is 3.7783e-08

## COMPUTER PROBLEMS 5.2

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.

(1) （（x^2*ln(x)在[1,3]上的积分）——公式复制不过来）   (2)（x/(sqrt(x^4)+1)在[0,1]上的积分）

The code of function Trapezoid Rule is:

%Composite Trapezoid Rule

function xc=ComTR(f,a,b,m)

h=(b-a)/m;

y0=f(a);ym=f(b);yk=0;

for i=1:m-1

x=a+i*h;

y=f(x);

yk=yk+y;

end

xc=h*(y0+ym+2*yk)/2;

the code of the problem is:

%5-2-1

clear;clc;

y1=@(x)x^2*log(x);a1=1;b1=3;

y2=@(x)x/sqrt(x^4+1);a2=0;b2=1;

x1_16=ComTR(y1,a1,b1,16);x1_32=ComTR(y1,a1,b1,32);

x2_16=ComTR(y2,a2,b2,16);x2_32=ComTR(y2,a2,b2,32);

the solution is: (1) we can have the exact value of the solution is 6.998621709124097

when n=16, the integral is: 7.00980923592,the error is about 0.0112; when n=32, the integral is 7.00141850616,the error is about 0.0028.

(2) when n=16, the integral is: 0.44036118262; when n=32,the integral is 0.44060540767.

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

Solution:

The code of the function Composite Simpson's Rule is:

%Composite Simpson's Rule

function xc=ComSR(f,a,b,m)

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

y0=f(a);y2m=f(b);

y1=0;y2=0;

for i=1:(2*m-1)

x(i)=a+i*h;

y(i)=f(x(i));

if mod(i,2)==1

y1=y1+y(i);

else

y2=y2+y(i);

end

end

xc=h*(y0+y2m+4*y1+2*y2)/3;

the code solving the question is:

%CP5-2-2

y1=@(x)x^2*log(x);a1=1;b1=3;

y2=@(x)x/sqrt(x^4+1);a2=0;b2=1;

x1_16=ComSR(y1,a1,b1,16);x1_32=ComSR(y1,a1,b1,32);

x2_16=ComSR(y2,a2,b2,16);x2_32=ComSR(y2,a2,b2,32);

the solution is: (1) when n=16, the integral is: 6.99862159624; when n=32, the integral is 6.99862170206

(2) when n=16, the integral is: 0.44068681602; when n=32,the integral is 0.44068679491.

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

(a)  exp(x^2)在[0,1]上的积分  (f)（cos(exp(x))在[0,pi]上的积分）

Solution:

The code of the problem is:

%5-2-3

clear;clc;

y1=@(x)exp(x^2);a1=0;b1=1;

y2=@(x)cos(exp(x));a2=0;b2=pi;

x1_16=ComTR(y1,a1,b1,16);x1_32=ComTR(y1,a1,b1,32);

x2_16=ComTR(y2,a2,b2,16);x2_32=ComTR(y2,a2,b2,32);

the solution is: (a) when n=16, the integral is: 1.464420310149482; when n=32, the integral is 1.463094102606428

(f) when n=16, the integral is: -0.278012950965426; when n=32,the integral is -0.356789537332442

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

Solution:

The code of the problem is:

%5-2-4

clear;clc;

y1=@(x)exp(x^2);a1=0;b1=1;

y2=@(x)cos(exp(x));a2=0;b2=pi;

x1_16=ComSR(y1,a1,b1,16);x1_32=ComSR(y1,a1,b1,32);

x2_16=ComSR(y2,a2,b2,16);x2_32=ComSR(y2,a2,b2,32);

the solution is: (a) when n=16, the integral is: 1.462652033425411; when n=32, the integral is 1.462651763901493

(f) when n=16, the integral is: -0.383048399454781; when n=32,the integral is -0.376327835672488.

## EX5.3

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

Solution:

From the steps of Romberg Integration we can have: h1=1, h2=1/2 and h3=1/4;

R11=(1+e)/2; R21=(1+e)/4+exp(1/2)/2;R31=(1+e)/8+(exp(1/2))/4+(exp(1/4)+exp(3/4))/4;

R22=((1+e)/2+2*exp(1/2))/3;R32=((1+e)/4+exp(1/2)/2+exp(1/4)+exp(3/4))/3;

R33=(7*(1+e)/6+2*exp(1/2)+16*(exp(1/4)+exp(3/4))/3)/15;

We also can solve the problem by programing the following code.

%Romberg Integration

function r=romberg(f,a,b,n)

h=(b-a)./(2.^(0:n-1));

r(1,1)=(b-a)*(f(a)+f(b))/2;

for j=2:n

t=0;

for i=1:2^(j-2)

t=t+f(a+(2*i-1)*h(j));

end

r(j,1)=r(j-1,1)/2+h(j)*t;

for k=2:j

r(j,k)=(4^(k-1)*r(j,k-1)-r(j-1,k-1))/(4^(k-1)-1);

end

end

the solution is r=[1.8591,0,0;1.7539,1.7189,0;1.7272,1.7183,1.7183]

# CHAPTER SIX: Ordinary Differential Equations

## EXERCISE6.1

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’=y*t^2  (c) y’=2(t+1)y  (d) y’=5*y*t^4  (e)y’=1/y^2  (f) y’=t^3/y^2

Solution:

By using the separation of variables we can have:

(a) y=(t^2)/2+c and with y(0)=1 we can have that c=1; then the solution of (a) is y=(t^2)/2+1;

(b) |y|=exp((t^2)/3+c) and with y(0)=1 we can have that c=0; then the solution of (b) is |y|= exp((t^2)/3)

(c) |y|=exp(t^2+2*t+c) and with y(0)=1 we can have that c=0; the solution of (c) is |y}=exp(t^2+2*t)

(d)|y|=exp(t^5+c) and with y(0)=1 we can have that c=0; it means the solution is |y|= exp(t^5)

(e)y^3=3*t+c and we can get c=1 with y(0)=1; the solution of (e) is y^3=3*t+1.

(f)y^3=3*t^4/4+c and we can get c=1 with y(0)=1; the solution of (f) is y^3=3*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.

Solution: when h=1/4, we can have that the values we choose on the interval [0,1] is t10=0; t1=1/4; t2=1/2; t3=3/4. Apply Euler’s Method to each equation we can have:

(a) y(0)=1, y’(0)=0; y(1/4)=y(0)+0*(1/4)=1, y’(1/4)=1/4; y(1/2)=y(1/4)+1/4*(1/4)=17/16, y’(1/2)=1/2; y(3/4)=y(1/2)+1/2*h=19/16, y’(3/4)=3/4; y(1)=y(3/4)+3/4*h=11/8.

We can get: w0=1; w1=1;w2=17/16;w3=19/16 and w4=11/8.

And we can have the exact value of y(1) is 3/2 from above problem 3. , then the error is 1/8=0.125

(b) y(0)=1, y’(0)=0; y(1/4)=y(0)+0*(1/4)=1, y’(1/4)=1/16; y(1/2)=y(1/4)+1/4*(1/4)=65/64, y’(1/2)=65/256; y(3/4)=y(1/2)+y’(1/2)*h=1.076171875, y’(3/4)= 0.6053466796875; y(1)=y(3/4)+y’(3/4)*h=1.227508544921875. And we can have the exact value of y(1) is exp(1/3)

We can get : w0=1,w1=1,w2=65/64,w3=1.076171875 and w4=1.227508544921875.

from above problem 3. ,then the error is 0.1681

(c)  y(0)=1, y’(0)=2; y(1/4)=y(0)+2*(1/4)=3/2, y’(1/4)=15/4; y(1/2)=39/16, y’(1/2)=117/16; y(3/4)= 273/64, y’(3/4)=1911/128; y(1)=y(3/4)+y’(3/4)*h=7.9980.

We can get that w0=1,w1=1,w2=13/8,w3=91/32 and w4=7.9980

And we can have the exact value of y(1) is  exp(3)  from above problem 3. , then the error is 12.0875

(d) y(0)=1, y’(0)=0; y(1/4)=y(0)+0*(1/4)=1, y’(1/4)=5/(4^4); y(1/2)=1+5/(4^5)= 1.0048828125, y’(1/2)=5*(1+5/(4^5))*(1/2)^4; y(3/4)=y(1/2)+y’(1/2)*h=1.0833, y’(3/4)= 1.7140; y(1)=1.5119. And we can have the exact value of y(1) is exp(1) from above problem 3. ,then the error is 1.2064

(e) y(0)=1, y’(0)=0; y(1/4)=y(0)+0*(1/4)=1, y’(1/4)=1; y(1/2)=y(1/4)+1/4*(1/4)=5/4, y’(1/2)=16/25; y(3/4)=y(1/2)+1/2*h=1.41, y’(3/4)=(100*100)/(141*141); y(1)=1.5357 . And from above problem 3. , then the error is 0.0517

(f) y(0)=1, y’(0)=0; y(1/4)= 1, y’(1/4)=1/64; y(1/2)=257/256, y’(1/2)=(256^2)/(8*257^2); y(3/4)=1.0310, y’(3/4)=0.3969; y(1)=y(3/4)+3/4*h=1.1302. And we can have the exact value of y(1) is 3/2 from above problem 3. ,then the error is 1/8=0.0748

## EXERCISE6.2

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   (d) y’=5*y*t^4  (e)y’=1/y^2  (f) y’=t^3/y^2

Solution:

(a)

w0=y0=1, y’(0)=0; w1=(1+1/32)=33/32, y’(1/4)=1/4;w2=33/32+3/32=9/8,y’(1/2)=1/2;w3=41/32, y’(3/4)=3/4;w4=23/16. And we have known that the exact value of y(1) is 3/2. Then we can have that the error is 1/16.

Then we use Matlab and the formula to solve the following equationsb-f:

%EX6.2.1

clear;

t=[0,1/4,1/2,3/4,1];h=1/4;

wb(1)=1;

for i=2:5

dy=t(i-1)*t(i-1)*wb(i-1);

dy1=t(i)*t(i)*(wb(i-1)+h*dy);

wb(i)=wb(i-1)+h*(dy+dy1)/2;

end

eb=abs(wb(5)-exp(1/3));

wc(1)=1;

for i=2:5

dy=2*(t(i-1)+1)*wc(i-1);

dy1=2*(t(i)+1)*(wc(i-1)+h*dy);

wc(i)=wc(i-1)+h*(dy+dy1)/2;

end

ec=abs(wc(5)-exp(3));

wd(1)=1;

for i=2:5

dy=5*(t(i-1)^4)*wd(i-1);

dy1=5*(t(i)^4)*(wd(i-1)+h*dy);

wd(i)=wd(i-1)+h*(dy+dy1)/2;

end

we(1)=1;ed=abs(wd(5)-exp(1));

for i=2:5

dy=1/(we(i-1)*we(i-1));

dy1=1/((we(i-1)+h*dy)*(we(i-1)+h*dy));

we(i)=we(i-1)+h*(dy+dy1)/2;

end

ee=abs(we(5)-nthroot(4,3));wf(1)=1;

for i=2:5

dy=t(i-1)^3/(wf(i-1)*wf(i-1));

dy1=t(i)^3/((wf(i-1)+h*dy)*(wf(i-1)+h*dy));

wf(i)=wf(i-1)+h*(dy+dy1)/2;

end

ef=abs(wf(5)-nthroot(7/4,3));

We can get that

wb= [1,1.0078,1.0477,1.1587,1.4054];   and the error eb=0.0097

wc= [1,1.7188,3.3032,7.0710,16.7935];  and the error ec= 3.2920

wd= [1,1.0024,1.0442,1.3077,2.7068];  and the error ed=0.0115

we= [1,1.2050,1.3570,1.4810,1.5871];  and the error ee=2.9956e-4

wf= [1,1.0020,1.0193,1.0823,1.2182];   and the error ef=0.0132