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

### 17074198 作业九

5.1Exercises

2. Problem: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:

(a) The three-point centered-difference formula evaluates to

f(0)(f(x+h)-f(x-h))/2h=(f(0.1)-f(-0.1))/0.2=5(e^0.1-1/e^0.1).

(b) The three-point centered-difference formula evaluates to

f(0)(f(x+h)-f(x-h))/2h=(f(0.01)-f(-0.01))/0.02=50(e^0.01-1/e^0.01)

(c) The three-point centered-difference formula evaluates to

f(0)(f(x+h)-f(x-h))/2h=(f(0.001)-f(-0.001))/0.002=500(e^0.001-1/e^0.001)

5.1Computer Problems

1. Problem:Make a table of the error of the three-point centered-difference formula for f(0), where f (x) = sin x cos x, 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?

Code:

function [df,err]=TPCDF(df1,f,h,x)

F=f(x+h)-f(x-h);

df=F/2*h;

err=df-df1(x);

Result:

>> syms x

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

>> df1=diff(f)

df1 =

cos(x) + sin(x)

>> df1=@(x) cos(x) + sin(x);

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

>> x=0;

>> h=10^-1;

>> [df,err]=TPCDF(df1,f,h,x)

df =

0.9983

err =

-0.0017

5.2 Exercises

1.Problem: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)(π/2,0)cosx dx

Solution:(b)

(π/2,0)cosx dx=sinx|(π/2,0)=1

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

(π/2,0)cosx dx(π/2/2)(y0+y1)

=(π/4)(cos0+cosπ/2)

=π/40.7854

err=1-0.7854=0.2146

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

(π/2,0)cosx dx(h/2)(y0+ym+2y1)

=(π/8)(cos0+cosπ/2+2cosπ/4)

=(π/8)(1+0+2*2/2)

0.9481

err=1-0.9481=0.0519

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

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

=(π/16)(cos0+cosπ/2+2(cosπ/8+cosπ/4+cos3π/8))

=(π/16)(1+0+2(0.9239+2/2+0.3827))

0.9871

err=1-0.9871=0.0129

2.Problem:Apply the Composite Midpoint Rule with m = 1,2, and 4 panels to approximate the integrals in Exercise 1, and report the errors.(b)(π/2,0)cosx dx

Solution:(b)

When m=1,h=b-a=π/2,

(π/2,0)cosx dx(π/2)f(w1)

=(π/2)*(2/2)

=2π/41.1107

err=1.1107-1=0.1107

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

(π/2,0)cosx dx(π/4)(f(w1)+f(w2))

=(π/4)(cos(π/8)+cos(3π/8))

1.0262

err=1.0262-1=0.0262

When m=4,h=π/8,

(π/2,0)cosx dx(π/8)(f(w1)+f(w2)+f(w3)+f(w4))

=(π/8)(cos(π/16)+cos(3π/16)+cos(5π/16)+cos(7π/16))

1.0065

err=1.0065-1=0.0065

3.Problem:Apply the composite Simpsons Rule with m = 1,2, and 4 panels to the integrals in Exercise 1,and report the errors.(b)(π/2,0)cosx dx

Solution:(b)

When m=1,h=π/4,

(π/2,0)cosx dx(π/4/3)(y0+y2+4y1)

=(π/12)(cos0+cos(π/2)+4cos(π/4))

1.0023

err=1.0023-1=0.0023

When m=2,h=π/8,

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

=(π/24)(cos0+cos(π/2)+4(cos(π/8)+cos(3π/8))+2cos(π/4))

1.0001

err=1.0001-1=0.0001

When m=4,h=π/16,

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

=(π/48)(cos0+cos(π/2)+4(cos(π/16)+cos(3π/16)+cos(5π/16)+cos(7/16))

+2(cos(π/8)+cos(π/4)+cos(3π/8))

1.0000

err=0

4.Problem:Apply the composite Simpsons Rule with m = 1,2, and 4 panels to the integrals, and report the errors.(b)(1,0)(1/(1+x^2))dx

Solution:(b)

(1,0)(1/(1+x^2))dx=arctanx|(1,0)0.7854

When m=1,h=1/2,

(1,0)(1/(1+x^2))dx(1/2/3)(y0+y2+4y1)

=(1/6)(1/(1+0^2)+1/(1+1)+4*(1/(1+1/4)))

=0.7833

err=0.7854-0.7833=0.0021

When m=2,h=1/4,

(1,0)(1/(1+x^2))dx(1/4/3)(y0+y4+4(y1+y3)+2y2)

=(1/12)(1/(1+0)+1/(1+1)+4(1/(1+(1/4)^2)+1/(1+(3/8)^2))+2(1/(1+(1/2)^2)))

0.8643

err=0.8643-0.7854=0.0789

When m=4,h=1/8,

(1,0)(1/(1+x^2))dx(1/8/3)(y0+y8+4(y1+y3+y5+y7)+2(y2+y4+y6))

=(1/24)(1/(1+0)+1/(1+1)+4(1/(1+1/64)+1/(1+9/64)+1/(1+25/64)+1/(1+49/64))

+2(1/(1+1/16)+1/(1+1/4)+1/(1+9/16)))

0.7854

err=0

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) (3,1)(x^2lnx)dx   (h)(1,0)(x/x^4+1)dx

Code:

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

h=(b-a)/m;

y0=f(a);

ym=f(b);

mult=0;

for i=1:m-1

yi=f(a+i*h);

mult=mult+yi;

end

y=(h/2)*(y0+ym+2*mult);

Result:

(d)

(3,1)(x^2lnx)dx =9ln(3)-26/9=6.9986

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

>> a=1;

>> b=3;

>> m=16;

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

y =

7.0098

err=7.0098-6.9986=0.0112

>> m=32;

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

y =

7.0014

err=7.0014-6.9986=0.0028

(h)

(1,0)(x/x^4+1)dx

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

>> a=0;

>> b=1;

>> m=16;

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

y =

0.4404

>> m=32;

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

y =

0.4406

2.Problem:Apply the composite Simpsons Rule to the integrals in Computer Problem 1. Use m = 16 and32, and report errors.(d)(3,1)(x^2lnx)dx   (h)(1,0)(x/x^4+1)dx

Code:

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

% f被积函数；a积分下限；b积分上限；m子区间个数（将x分为多少个区间）

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

s1=0; s2=0;

for i=1:m

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

s1=s1+f(x);

end

for j=1:(m-1)

x=a+2*j*h;

s2=s2+f(x);

end

y=h/3*(f(a)+2*s2+4*s1+f(b));

Result:

(a)

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

>> a=1;

>> b=3;

>> m=16;

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

y =

6.9986

>> m=32;

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

y =

6.9986

(b)

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

>> a=0;

>> b=1;

>> m=16;

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

y =

0.4407

>> m=32;

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

y =

0.4407

3.Problem:Use the composite Trapezoid Rule with m = 16 and 32 panels to approximate the definite integral.(a)(1,0)e^x^2dx   (f)(π,0)cos(e^x)dx

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

Result:

(a)

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

trapezoid(0,1,16,f)

ans =

1.464420310149482

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

ans =

1.463094102606428

(f)

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

trapezoid(0,pi,16,f)

ans =

-0.278012950965426

ans =

-0.2780

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

ans =

-0.356789537332442

4.Problem:Apply the composite Simpsons Rule to the integrals of Computer Problem 3, using m = 16 and 32.(a)(1,0)e^x^2dx   (f)(π,0)cos(e^x)dx

Code:

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

% f被积函数；a积分下限；b积分上限；m子区间个数（将x分为多少个区间）

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

s1=0; s2=0;

for i=1:m

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

s1=s1+f(x);

end

for j=1:(m-1)

x=a+2*j*h;

s2=s2+f(x);

end

y=h/3*(f(a)+2*s2+4*s1+f(b));

(a)

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

y=Simpson1(f,0,1,16)

y =

1.462652033425411

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

y =

1.462651763901493

(f)

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

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

y =

-0.123445767930410

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

y =

-0.123445911966419

5.3Exercises

1. Problem:Apply Romberg Integration to find R33 for the integrals.

(c)(1,0)e^xdx

Solution:(c)

R11=(h1/2)(f(a)+f(b))=(1+e)/2

R21=R11/2+h2f((a+b)/2)=(1+e)/4+(e^(1/2))/2

R31=R21/2+h3(f(a+h3)+f(a+h3))=(1+e)/8+(e^(1/2)+(e^(1/4)+e^(3/4))/4

R32=(3R31-R21)/3=(1/3)((1+e)/4+(e^(1/2))/2+e^(1/4)+e^(3/4))

R22=(4R21-R11)/3=(1+e)/6+2(e^(1/2))/3

R33=(1/15)(16R32-R22)=4(1+e)/45-(1+e)/6-22(e^(1/2))/45+16(e^(1/4)+e^(3/4))/45=-0.0777

6.1Execises

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

Solution:

(a) dy/dt=t

∫dy=tdt

y=t^2/2+k

y(0)=1=k

y=t^2/2+1

(b) dy/dt=t^2*y

dy/y=t^2dt

ln|y|=t^3/3+k

|y|=e^(t^3/3+k)

y(0)=1=e^k

y=e^(t^3/3)

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

dy/y=2(t+1)dt

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

|y|=e^(t^2+2t+k)

y(0)=1=e^k

y=e^(t^2+2t)

(d) dy/dt=5t^4*y

dy/y=5t^4dt

ln|y|=t^5+k

|y|=e^(t^5+k

y(0)=1=e^k

y=e^(t^5)

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

y^2dy=dt

y^3/3=t+k

y=(3t+k)^(1/3)

y(0)=1=k^(1/3)

k=1

y=(3t+1)^(1/3)

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

y^2dy=t^3dt

y^3/3=t^4/4+k

y=(3t^4/4+k)^(1/3)

y(0)=1=k^(1/3)

k=1

y=(3t^4/4+1)^(1/3)

5.Problem:Apply Eulers 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 find the error at t = 1 by comparing with the correct

solution.

(a)y=t   (b)y=t^2*y   (c)y=2(t+1)y   (d)y=5t^4y   (e)y=1/y^2  (f)y=t^3/y^2

Solution:

(a) w0=1

w1=w0+ht=1+1/4*1/4=17/16

w2=w1+ht=17/16+1/4*1/2=19/16

w3=w2+ht=19/16+1/4*3/4=11/8

w4=w3+ht=11/8+1/4*1=13/8

y(1)=3/2,err=w4-y(1)=0.1250

(b) w0=1

w1=w0+h*t^2*w0=1+(1/4)*(1/16)=65/64

w2=w1+h^t^2*w1=65/64+(1/4)*(1/4)*(65/64)=1.0791

w3=w2+h^t^2*w2=1.0791+(1/4)*(9/16)*1.0791=1.2308

w4=w3+h*t^2*w3=1.2308+(1/4)*1*1/2308=1.2309

y(1)=e^(1/3),err=y(1)-w4=e^(1/3)-1.2309

(c) w0=1

w1=w0+h*2(t+1)*w0=1+(1/4)*2*(5/4)*1=13/8

w2=w1+h*2(t+1)*w1=13/8+(1/4)*2*(3/2)*(13/8)=2.8438

w3=w2+h*2(t+1)*w2=2.8438+(1/4)*2*(7/4)*2.8438=5.3321

w4=w3+h*2(t+1)*w3=5.3321+(1/4)*2*2*5.3321=10.6642

y(1)=e^3,err=e^3-10.6642

(d) w0=1

w1=w0+h*5t^4*w0=1+(1/4)*5*(1/4)^4*1=1.0049

w2=w1+h*5t^4*w1=1.0049+(1/4)*5*(1/2)^4*1.0049=1.0834

w3=w2+h*5t^4*w2=1.0834+(1/4)*5*(3/4)^4*1.0834=1.5119

w4=w3+h*5t^4*w3=1.5119+(1/4)*5*(1)^4*1.5119=3.4018

y(1)=e,err=e-3.4018

(e) w0=1

w1=w0+h*(1/w0^2)=1+(1/4)*1=5/4

w2=w1+h*(1/w1^2)=5/4+(1/4)*(16/25)=1.4100

w3=w2+h/(w2^2)=141/100+(1/4)*(10000/141*141)=1.5357

w4=w3+h/(w3^2)=1.6417

err=|y4-w4|=|4^(1/3)-1.6417|=0.0543

(f) w0=1

w1=w0+h*(t0^3)*(1/w0^2)=1+0=1
w2=w1+h*(t1^3)*(1/w1^2)=1+(1/4)*(1/64)=257/256=1.0039
w3=w2+h*(t2^3)*(1/w2^2)=257/256+(1/4)*(8/64)*(256*256/257*257)=1.0349
w4=w3+h*(t3^3)*(1/w3^2)=1.1334
e4=|y4-w4|=|(7/4)^(1/3)-1.1334|=0.0717

6.2Exercises

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

Solution:

(a) w0=1

w1=w0+(h/2)*(2t+h)=1+(1/8)*(1/4)=1.0313

w2=w1+(h/2)*(2t+h)=1.0313+(1/8)*(2*(1/4)+(1/4))=1.1251

w3=w2+(h/2)*(2t+h)=1.1251+(1/8)*(2*(1/2)+(1/4))=1.2813

w4=w3+(h/2)*(2t3+h)=1.2813+(1/8)*(2*(3/4)+(1/4))=1.5001

err=0.0001

(b) w0=1
w1=w0+(h/2)*((t0^2)*y0+((t0+h)^2)*(w0+h*(t0^2)*y0))=1.0078
w2=w1+(h/2)*((t1^2)*y1+((t1+h)^2)*(w1+h*(t1^2)*y1))=1.0477
w3=w2+(h/2)*((t2^2)*y2+((t2+h)^2)*(w2+h*(t2^2)*y2))=1.1587
w4=w3+(h/2)*((t3^2)*y3+((t3+h)^2)*(w3+h*(t3^2)*y3))=1.4054
e4=|y4-w4|=0.0097

(c) w0=1
w1=w0+(h/2)*(2y0*(t0+1)+2(t0+h+1)*(w0+2h(t0+1)y0))=1.7188
w2=w1+(h/2)*(2y1*(t1+1)+2(t1+h+1)*(w1+2h(t1+1)y1))=3.3032
w3=w2+(h/2)*(2y2*(t2+1)+2(t2+h+1)*(w2+2h(t2+1)y2))=7.0710
w4=w3+(h/2)*(2y3*(t3+1)+2(t3+h+1)*(w3+2h(t3+1)y3))=16.7935
e4=|y4-w4|=3.2920

(d) w0=1

w1=w0+(h/2)*(5*t0^4*w0+5*((t0+(1/4)^4)*(w0+(1/4)*(5t0^4*w0))

=1+(1/8)*(0+5*(1/4)^4*1)=1.0024

w2=w1+(1/8)*(5*(1/4)^4*w1+5*(1/2)^4*(w1+(1/4)*5*(1/4)^4*w1))=1.0442

w3=w2+(1/8)*(5*(1/2)^4*w2+5*(3/4)^4*(w2+(1/4)*5*(1/2)^4*w2))=1.2575

w4=w3+(1/8)*(5*(3/4)^4*w3+5*(1)^4*(w3+(1/4)*5*(3/4)^4*w3))=2.6030

err=|2.6030-e|=0.1153

(e) w0=1

w1= w0+h/2(f(t0,w0)+f(t0+h,w0+hf(t0,w0)))=1.205000000000000

w2= w1+h/2(f(t1,w1)+f(t1+h,w1+hf(t1,w1)))=1.356993862239376

w3= w2+h/2(f(t2,w2)+f(t2+h,w2+hf(t2,w2)))=1.480971731811849

w4= w3+h/2(f(t3,w3)+f(t3+h,w3+hf(t3,w3)))=1.587101493145249

e=|y-w4|=2.995588229504076e-04

(f) w0=1

w1= w0+h/2(f(t0,w0)+f(t0+h,w0+hf(t0,w0)))=1.001953125000000

w2= w1+h/2(f(t1,w1)+f(t1+h,w1+hf(t1,w1)))=1.019342601468375

w3= w2+h/2(f(t2,w2)+f(t2+h,w2+hf(t2,w2)))=1.082264953617726

w4= w3+h/2(f(t3,w3)+f(t3+h,w3+hf(t3,w3)))=1.218241931698513

e=|y-w4|=0.013170799610898