Homework 4.4: Numerical Calculus

Andy Long, Spring, 2024

pp. 172, 1acf, 2acf, 3acf, 7, 8, 11ab, 12, 15, 27, 28, 29

#1acf

compMidpoint[f_,a_,b_,n_]:=
Module[{h=(b-a)/n,int},
int=0;
For[i=1/2,i<=n,i++,
int=int+f[a+i h];
];
int*h
]

1a.

f[x_]:=Log[Sin[x]]
{a,b}={1.0,3.0};
compMidpoint[f,a,b,3]
Integrate[f[x],{x,a,b}]

1c.

f[x_]:=E^x Log[x]/x
{a,b}={1.0,4.0};
compMidpoint[f,a,b,3]
Integrate[f[x],{x,a,b}]

1f.

f[x_]:=(x^2-1)/(x^2+1)
{a,b}={0.0,1.0};
compMidpoint[f,a,b,3]
Integrate[f[x],{x,a,b}]

#2acf

compTrap[f_,a_,b_,n_]:=
Module[{h=(b-a)/n,fa=f[a],fb=f[b],int},
int=(fa+fb)/2;
For[i=1,i<=n-1,i++,
int=int+f[a+i h];
];
int*h
]

1a.

f[x_]:=Log[Sin[x]]
{a,b}={1.0,3.0};
compTrap[f,a,b,3]
Integrate[f[x],{x,a,b}]

1c.

f[x_]:=E^x Log[x]/x
{a,b}={1.0,4.0};
compTrap[f,a,b,3]
Integrate[f[x],{x,a,b}]

1f.

f[x_]:=(x^2-1)/(x^2+1)
{a,b}={0.0,1.0};
compTrap[f,a,b,3]
Integrate[f[x],{x,a,b}]

#3acf

Note that we can once again get off cheaply, by using two approximations, and averaging them together in a sensible way:

compSimp[f_,a_,b_,n_]:=(2 compMidpoint[f,a,b,n]+compTrap[f,a,b,n])/3

1a.

f[x_]:=Log[Sin[x]]
{a,b}={1.0,3.0};
compSimp[f,a,b,3]
Integrate[f[x],{x,a,b}]

1c.

f[x_]:=E^x Log[x]/x
{a,b}={1.0,4.0};
compSimp[f,a,b,3]
Integrate[f[x],{x,a,b}]

1f.

f[x_]:=(x^2-1)/(x^2+1)
{a,b}={0.0,1.0};
compSimp[f,a,b,3]
Integrate[f[x],{x,a,b}]

#7 & 8

#7: 109 will suffice (turns out 3 is enough) (Check -- with p. 294)

f[x_]:=Sin[x]^4
{a,b}={0.0,Pi};
true=Integrate[f[x],{x,0,Pi}]
estimate=compTrap[f,a,b,1]
Manipulate[
{compTrap[f,a,b,n],compMidpoint[f,a,b,n],true*1.0},
{n,1,20,1}
]
Solve[Abs[true-estimate]/n^2==10^(-4),n]

#8: 141 will suffice (turns out 3 is enough)

estimate=compMidpoint[f,a,b,1]
Solve[Abs[true-estimate]/n^2==10^(-4),n]

#11ab

#12

f[x_]:=1/(1-x)
f4[x_]=D[f[x],{x,4}]
Plot[Abs[f4[x]],{x,2,4},PlotRange->All]
Mbound=f4[2]
h=0.1;
Abs[(4-2)/90 h^4 Mbound]

#15

I’m starting with adaptive trapezoidal, because I’ve already got it in hand. I check it first against the author’s example; then adapt it for Simpson’s.

aTraptiveQuad[f_,enda_,endb_,fenda_,fendb_,eps_,maxits_,verbose_]:=
Module[{a=enda,b=endb,fa=fenda,fb=fendb,h,c,fc,err,int},
h=(b-a)/2;
c=a+h;
fc=N[f[c]]; (* 
(T2-T1)/3 
= 
h((fa+2fc +fb)/2-(fa+fb))/3 
= 
-h/6(fa-2fc+fb)*)
err=h /6 *(fa -2 fc+ fb); 
If[verbose,Print[{a,b,h(fa+2fc +fb)/2,h(fa+fb),Abs[err]}]];
If[Abs[err]<eps || maxits == 0
,int=h(fa+4fc +fb)/3;(* Simpson's! Ironic, no? *)
,
int=
aTraptiveQuad[f,a,c,fa,fc,eps/2,maxits-1,verbose]
+
aTraptiveQuad[f,c,b,fc,fb,eps/2,maxits-1,verbose]
];
Return[int]
]

Check that it’s working against the table of page 172.

Note that my approximation is a little better than the author’s, because I compute Simpson’s as my estimate rather than T2 -- since I’ve got everything I need to produce it.

f[x_]:=Log[x+3]
{a,b}={0.0,3.0};
eps=0.006;
its=30;
NIntegrate[f[x],{x,a,b}]
aTraptiveQuad[f,a,b,f[a],f[b],eps,its,True]

aSimptiveQuad[f_,enda_,endb_,fenda_,fmid_,fendb_,eps_,maxits_,verbose_]:=
Module[{
a=enda,b=endb,c=Mean[{enda,endb}],
fa=fenda,fm=fmid,fb=fendb,
fml,fmr,
ends,mids,h,halfh,err,int},

h=(b-a)/2; (* NOTE: h is the half width of the interval. *)
halfh=h/2; (* NOTE: halfh is the half half-width of the interval!:) *)
fml=f[a+halfh];
fmr=f[b-halfh]; 
(*
s1=h/3(fa+4fm+fb);
s2=h/6(fa+4(fml+fmr)+2fm+fb); 
err=(s2-s1)/15;
*)
(* These economize the calculations above: *)
ends=fa+fb;
mids=4(fml+fmr);
err=(ends-mids+6fm)*h/90; 
If[verbose, (* I have to actually do the calculations of s1 and s2: *)Print[{a,b,h/3(fa+4fm+fb),h/6(fa+4(fml+fmr)+2fm+fb),Abs[err]}]];

If[Abs[err]<eps || maxits == 0, (* We've hit the end of the line: *)
int=(ends+mids+2fm)h/6, (* This is s2 *)
int= (* Otherwise split divide the error in half, and split the interval: *)
aSimptiveQuad[f,a,c,fa,fml,fm,eps/2,maxits-1,verbose]
+ (* not forgetting the reduce the number of iterations possible... *)
aSimptiveQuad[f,c,b,fm,fmr,fb,eps/2,maxits-1,verbose]
];

Return[int]
]

Our author means for you to go through the calculations, I guess, as one can see from the solution in the back of the book:

f[x_]:=Log[Sin[x]]
{a,b}={1.0,3.0};
eps=.002;
its=30;
NIntegrate[f[x],{x,a,b}]
aSimptiveQuad[f,a,b,f[a],f[(a+b)/2.0],f[b],eps,its,True]

#27

#28

f[x_]:=Log[Sin[x]]
{a,b}={1.0,3.0};
eps=.002;
its=30;
NIntegrate[f[x],{x,a,b}]
aTraptiveQuad[f,a,b,f[a],f[b],eps,its,True]
aSimptiveQuad[f,a,b,f[a],f[(a+b)/2.0],f[b],eps,its,True]

#29

f[x_]:=Log[x+1]
{a,b}={0.0,1.0};
eps=10^(-4.0);
its=30;
NIntegrate[f[x],{x,a,b}]
aTraptiveQuad[f,a,b,f[a],f[b],eps,its,False]
aSimptiveQuad[f,a,b,f[a],f[(a+b)/2.0],f[b],eps,its,False]