﻿ 17074199 5.1+5.2+5.3+6.1+6.2 作业九 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 17074199 5.1+5.2+5.3+6.1+6.2 作业九

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.

(a) According to the Taylor method,the three-point centered-difference formula evaluates to

f’(x)=[f(x+h)-f(x-h)]/2h-(h^2/6)*f’’’(c)

=[e^(0.1)-e^(-0.1)]/(2*0.1)-[(0.1)^2/6]*f’’’(c)   -0.1<c<0.1

The approximate solution is f’(x) [f(x+h)-f(x-h)]/2h

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

Therefore,the approximate solution is 1.007

(b) According to the Taylor method ,The three-point centered-difference formula evaluates to  f’(x)=[f(x+h)-f(x-h)]/2h-(h^2/6)*f’’’(c) -0.01<c<0.01

The approximate solution is f’(x) [f(x+h)-f(x-h)]/2h

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

Therefore,the approximate solution is 1.0000167.

(c) According to the Taylor method ,The three-point centered-difference formula evaluates to f’(x)=[f(x+h)-f(x-h)]/2h-(h^2/6)*f’’’(c) -0.001<c<0.001

The approximate solution is f’(x) [f(x+h)-f(x-h)]/2h

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

Therefore,the approximate solution is 1.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) = sinxcosx, 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?

According to the Taylor method,the three-point centered-difference formula evaluates to

f’(x)=[f(x+h)-f(x-h)]/2h-(h^2/6)*f’’’(c)

The three-point formula (5.7) yields f’(x)[f(x+h)-f(x-h)]/2h[sin(x+h)-cos(x+h)-sin(x-h)+cos(x-h)]/2h

The results of these formulas for x=0 and a wide range of the increment size h,along with errors compared with the correct value sin0-cos0=-1,are given in the following table:

>> syms h

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

h =

0.100000000000000

0.010000000000000

0.001000000000000

0.000100000000000

0.000010000000000

0.000001000000000

0.000000100000000

0.000000010000000

0.000000001000000

0.000000000100000

0.000000000010000

0.000000000001000

>> y=[sin(h)-cos(h)-sin(-h)+cos(-h)]/(2*h)

y =

0.998334166468282

0.099998333341667                   0

0.009999998333333                   0

0.000999999998333                   0

0.000099999999998                   0

0.000010000000000                   0

0.000000999999999                   0

0.000000100000001                   0

0.000000010000000                   0

0.000000001000000                   0

0.000000000100000                   0

0.000000000010000                   0

Error=1-y

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)π/2

0     cosx dx

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

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

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

(cosx)(0,π/2)=(π/2)/2(y0+y1)= [(π/2)/2](y0+y1)= (π/4)(cos0+cosπ/2)=(π/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*cos(π/4))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,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. Apply the Composite Midpoint Rule with m = 1,2, and 4 panels to approximate the integrals

in Exercise 1, and report the errors.

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

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

The Composite Midpoint Rule:

(cosx)(0, π/2)=hf(w)+((b-a)*h^2/24)*f’’(c)     w=x0+h/2

(1) When m=1,h=(b-a)/m=(π/2-0)/1=b-a=π/2, the midpoints are w1=π/4.

(π/2,0)cosx dx=h∑f(xm) =hf(w1)(π/2)f(w1)

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

=2π/41.1107

err=1.1107-1=0.1107

(2) When m=2,h=(b-a)/m=(π/2-0)/2=π/4, the midpoints are w1=π/8,w2=3π/8.

(π/2,0)cosx dx=h∑f(xm)= h∑[f(w1)+f(w2)](π/4)(f(w1)+f(w2))

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

1.0262

error=1.0262-1=0.0262

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

the midpoints are w1=π/16,w2=3π/16,w3=5π/16,w4=7π/16.

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

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

1.0065

Error=1.0065-1=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.

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

the composite Simpson’s Rule is f(x)|b a=(h/3)[y0+y2m+4*y2i-1(i=1m)+2y2i(i=1m-1)-(b-a)h^4/180*f(4)(c)

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

(π/2,0)cosx dx(π/4/3)(y0+y2+4y1) .So the panels are y0=f(0),y1=f(π/4),y2=f(π/2).

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

error=|1-1.0022799|=0.0022799

(2)    When m=2, h=(b-a)/2m=(π/2-0)/4=π/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.0001346

error=|1-1.0001346|=0.0001346.

(3)    When m=4,h==(b-a)/2m=(π/2-0)/8=π/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.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)

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

(1) When m=1,h=(b-a)/2m=(1-0)/2=1/2, So the panels are y0=0,y1=f(1/2),y2=f(1).

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

=47/600.783333

err=0.7854-0.7833=0.0021

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

So the panels are y0=0,y1=f(1/4),y2=f(1/2),y3=f(3/4),y4=f(1).

(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

(3) When m=4,h==(b-a)/2m=(1-0)/8=1/8, So the panels are

y0=f(0),y1=f(1/8),y2=f(1/4),y3=f(3/8),y4=f(1/2),y5=f(5/8),y6=f(3/4),y7=f(7/8),y8=f(1).

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

The definite integral is 7.0098 and 7.0014.

(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)    (x^2)*(log(x));

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

ans =

7.0098

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

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

ans =

7.0014

The real solution of (x^2) (ln(x))(1,3)is 9ln3-26/9.

m=16,error=9ln3-26/9-7.0098;

m=32,error=9ln3-26/9-7.0014.

The definite integral is 0.4404 and 0.4406.

(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)  x/((x^4+1)^(1/2));

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

ans =

0.4404

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

ans =

0.4406

The real solution of x/((x^4+1)^(1/2))(0,1) is 1/2(ln(2^(1/2))+1)

m=16,error=1/2(ln(2^(1/2))+1)-0.4404

m=32, error=1/2(ln(2^(1/2))+1)-0.4406.

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

(d) (x^2) (ln(x))(1,3)    (h) x/((x^4+1)^(1/2))(0,1)

//错误

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

simpson=0;

for i=1:m

simpson1=simpson+f(a+(2*i-1)*(b-a)/(2*m));

end

for j=1:m-1

simpson2=simpson+f(a+(2*j)*(b-a)/(2*m));

end

simpson=((b-a)/3*2*m)*(f(a)+f(b)+4*simpson1+2*simpson2);

end

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

(d)The real solution of (x^2) (ln(x))(1,3)is 9ln3-26/9.

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

y =

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

y =

6.998621702062214

m=16,error=9ln3-26/9-6.998621596247110;

m=32,error=9ln3-26/9-6.998621702062214

(h) The real solution of x/((x^4+1)^(1/2))(0,1) is 1/2(ln(2^(1/2))+1)

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

f =

包含以下值的 function_handle:

@(x)x/((x^4+1)^(1/2))

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

y =

0.440686816029813

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

y =

m=16,error=1/2(ln(2^(1/2))+1)-0.440686816029813=

m=32, error=1/2(ln(2^(1/2))+1)-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.464420310149482

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

ans =

1.463094102606428

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

ans =

-0.2780

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

ans =

-0.356789537332442

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

and 32.

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

>> 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.3 Exercises

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

dy/dt=t,dy=tdt,y=1/2t^2+C;

y(0)=0+C=1;

C=1;

y=1/2*t^2+1.

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

In|y|=(1/3)*t^3+C |y|=e^[(1/3)*t^3+C]

y(0)=1  C=0

dy/dt=2(t+1)y

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

In|y|=(t+1)^2+C

y=e^[(t+1)^2+C]

y(0)=1 C=-1

So y=e^[(t+1)^2-1]

dy/dt=5(t^4)y,(1/y)dy=(5t^4)dt;

ln|y|=t^5+C;

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

y(0)=e^C=1,C=0;

So,y=e^(t^5).

(e)y’=1/y^2

dy/dt=1/y^2

y^2dy=dt

(1/3)*y^3=t+C  y=[3(t+C)]^(1/3)

y(0)=1  C=1/3

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

(f)y’=t^3/y^2

dy/dt=t^3/y^2

y^2dy=t^3dt

(1/3)*y^3=(1/4)*y^4+C y=[(3/4)*y^4+3C]^(1/3)

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

y=(3/4*t^4+1)^(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 w i ,i = 0,...,4, and find the error at t = 1 by comparing with the correct

solution.

Euler method:

w0=y0    w(i+1)=wi+hf(ti,wi)  wi is approximate solution

h=1/4

(a)      y’=t       y=(1/2)*t^2+1.

Y0=1 y1=33/32 y2=9/8 y3=41/32 y4=3/2

h=1/4   ti=0,1/4,1/2,3/4,1

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

w(i+1)=wi+hti

w0=y0=1

w1=w0+(1/4)(t0)=1

w2=w1+(1/4)(t1)=1+(1/4)(1/4)=17/16

w3=w2+(1/4)(t2)=17/16+(1/4)(1/2)=19/16

w4=w3+(1/4)(t3)=19/16+(1/4)(3/4)=22/16

w(i+1)=wi+h*ti
w0=y0=1
w1=w0+1/4*t0=1
w2=w1+1/4*t1=1+1/16=1.0625
w3=w2+1/4*t2=17/16+2/16=1.1875
w4=w3+1/4*t3=19/16+3/16=1.3750
e4=|y4-w4|=0.125

 步数 ti wi yi ei 0 0 1 1 0 1 1/4 1 33/32 2 1/2 17/16 9/8 3 3/4 19/16 41/32 4 1 22/16=1.375 3/2=1.5 0.125

(b)y’=t^2*y   y=e^[(1/3)*ti^3]

h=1/4   ti=0,1/4,1/2,3/4,1

Y0=1 y1=e^(1/192) y2=e^(1/24) y3=e^(9/64) y4=e^(1/3)

W(i+1)=wi+hf(ti,wi)=wi+h(ti)^2*wi=wi*[1+h*(ti)^2]

w0=y0=1

w1=w0*[1+(1/4)*0]=1

w2=w1[1+(1/4)*(1/16)]=65/64

w3=w2[1+(1/4)*(1/4)]=(65/64)*(17/16)=1105/1024= 1.0791015625

w4=w3[1+(1/4)*(9/16)]=(1105/1024)*(73/64)=80665/65536= 1.230850219726563

w(i+1)=wi+h*(ti^2)*wi
w0=1
w1=w0+h*0=1
w2=w1+h*(t1^2)*w1=65/64=1.0156
w3=w2+h*(t2^2)*w2=65/64+(1/4)*(4/16)*(65/64)=1.0791
w4=w3+h*(t3^2)*w3=1105/1024+(1/4)*(9/16)*(1105/1024)=1.2309
e4=|y4-w4|=|e^(1/3)-1.2309|=0.1647

 步数 ti yi wi ei 0 0 1 1 0 1 1/4 1 1 2 1/2 e^(1/24) 65/64 3 3/4 e^(9/24) 1.0791015625 4 1 e^(1/3) 1.230850219726563

(c)y’=2(t+1)*y  y=e^[(t+1)^2-1]

h=1/4   ti=0,1/4,1/2,3/4,1

y0=1 y1= e^(9/16) y2= e^(5/4) y3=e^(-7/16) y4=e^3

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

W(i+1)=wi+h[2(ti+1)*wi]=[1+2h(ti+1)]wi

W0=y0=1

W1=1+2*1/4*1=3/2

W2=(3/2)*[1+1/2*5/4]=39/16

W3=(39/16)*[1+1/2*3/2]=(39/16)*(7/4)=273/64

W4=(273/64)*[1+1/2*7/4]= (273/64)*(7/8)=1911/512=3.732421875

Therefore,w(i+1)=wi+2h(ti+1)wi
w0=1
w1=w0+2h(t0+1)w0=1+1/2*1=1.5
w2=w1+2h(t1+1)w1=3/2+(1/2)*(3/2)*(5/4)=39/16=2.4375
w3=w2+2h(t2+1)w2=39/16+(1/2)*(39/16)*(6/4)=273/64=4.2656
w4=w3+2h(t3+1)w3=273/64+(1/2)*(273/64)*(7/4)=7.9980
e4=|y4-w4|=|e^3-7.9980|=12.0875

 步数 ti yi wi ei 0 0 1 1 0 1 1/4 e^(9/16) 3/2 2 1/2 e^(5/4) 39/16 3 3/4 e^(-7/16) 273/64 4 1 e^(3) 3.732421875

(d)y’=5*t^4*y   y=e^(t^5).

h=1/4   ti=0,1/4,1/2,3/4,1

Y0=1 y1=e^(1/1024)  y2=e^(1/32) y3=e^(273/1024) y4=e

w(i+1)=wi+hf(ti,wi)=wi+5h*(ti^4)*wi=[1+5h(ti^4)]wi

W0=y0=1

W1=1+(5/4)*0=1

W2=[1+(5/4)*(1/256)]=1+5/1024=1029/1024

W3=(1029/1024)*[1+5/4*1/16]= (1029/1024)*[69/64]=

W4=

Therefore,w(i+1)=wi+5h(ti^4)wi
w0=1
w1=w0+5h(t0^4)w0=1
w2=w1+5h(t1^4)w1=1+5*(1/4)*(1/256)*1=1029/1024=1.0049
w3=w2+5h(t2^4)w2=1029/1024+(5/4)*(16/256)*(1029/1024)=1.0834
w4=w3+5h(t3^4)w3=71001/65536+(5/4)*(81/256)*(71001/65536)=1.5119
e4=|y4-w4|=|e-1.5119|=1.2064

(e) y’=1/y^2   y=[3(t+1/3)]^(1/3)

h=1/4   ti=0,1/4,1/2,3/4,1

Y0=1 y1=e^(1/1024)  y2=e^(1/32) y3=e^(273/1024) y4=e

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

w(i+1)=wi+h/(wi^2)
w0=1
w1=w0+h/(w0^2)=1+1/4=5/4=1.25
w2=w1+h/(w1^2)=5/4+(1/4)*(16/25)=141/100=1.41
w3=w2+h/(w2^2)=141/100+(1/4)*(10000/141*141)=1.5357
w4=w3+h/(w3^2)=1.6417
e4=|y4-w4|=|4^(1/3)-1.6417|=0.0543

(f) y’=t^3/y^2  y=(3/4*t^4+1)^(1/3)

wi+1=wi+hf(ti,wi)
w(i+1)=wi+h*(ti^3)*(1/wi^2)
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.2 Exercises

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

approximation w 0 ,...,w 4 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 ′ = 5*t^4*y  (e) y ′=1/y^2    (f) y ′= t^3 /y^2

h=1/4, t0=0,t1=1/4,t2=1/2,t3=3/4,t4=1
(a)
w(i+1)=wi+(h/2)*(2ti+h)
w0=1
w1=w0+(h/2)*(2t0+h)=1.0313
w2=w1+(h/2)*(2t1+h)=1.1250
w3=w2+(h/2)*(2t2+h)=1.2813
w4=w3+(h/2)*(2t3+h)=1.5
e4=|y4-w4|=0

(b)
y=e^(t^3/3)
y0=1,y1=e^(1/192),y2=e^(8/192),y3=e^(27/192),y4=e^(1/3)
w(i+1)=wi+(h/2)*((ti^2)*yi+((ti+h)^2)*(wi+h*(ti^2)*yi))
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)
y=e^(t^2+2t)
y0=1,y1=e^(1/32),y2=e^(1/4),y3=e^(12/32),y4=e^3
w(i+1)=wi+(h/2)*(2yi*(ti+1)+2(ti+h+1)*(wi+2h(ti+1)yi))
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)

y= y’=5*t^4*y   y=e^(t^5).

h=1/4   ti=0,1/4,1/2,3/4,1

y0=1 y1=e^(1/1024)  y2=e^(1/32) y3=e^(273/1024) y4=e

(d)

w0=1

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

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

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

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

e= |y-w4|=0.011488678936460

(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