Numerical Analysis: Sections 1.3-1.4
Andy Long

Section 1.3: Algorithms

One might define algorithms as procedures for solving problems. Here’s one from Babylon, circa 2000 BCE (p. 13):

c=4
f[x_]:=c/x

As an iteration:

x0=1
y0=f[x0]
x1=(x0+y0)/2
y1=f[x1]
N[Abs[x1-y1]]
x2=(x1+y1)/2
y2=f[x2]
N[Abs[x2-y2]]
x3=(x2+y2)/2
y3=f[x3]
N[Abs[x3-y3]]
x4=(x3+y3)/2
y4=f[x4]
N[Abs[x4-y4]]

As a recursion (a fixed point problem -- the solution is x=g(x)):

g[x_]:=(x+f[x])/2
g[g[g[g[x0]]]]

xmax=2.5*c;
Plot[
{f[x],x}
,{x,0,xmax}
, AxesOrigin->{0,0}
, AspectRatio->1
, PlotRange->{{0,xmax},{0,xmax}}
]

Computer problem #1, p. 15:

it:=( x = (x + c/x)/2)

x=8
it
it
it
it
it
N[%]

The authors encourage us to experiment with lists, and the “%” command:

x=8
{it,it,it,it,it,it}
N[%]
N[%%,10]

I don’t mind the lists, but the “%” command makes me a little nervous, because if we introduce an intervening command, the meaning of “%” changes -- so this change will propogate through future commands, without another user being able to detect the intentions of the original author. Hence I prefer creating a new variable, with a sensible name, to store the result:

x=8
AppoximateRoots={it,it,it,it,it,it}
N[AppoximateRoots]
N[AppoximateRoots,20]

In addition to what the authors suggest we do, why don’t we have a look at the convergence in the values?

xs=AppoximateRoots;
ys=c/AppoximateRoots; (* f[xs] *)
Data=Transpose[{xs,ys}]
xmin=Min[Min[xs],Min[ys]];
xmax=Max[Max[xs],Max[ys]];

Plot1=
Plot[
{f[x],x}
,{x,.9*xmin,1.1*xmax}
, AspectRatio->1
, PlotRange->{{.9*xmin,1.1*xmax},{xmin,xmax}}
];
Plot2=
ListPlot[
Data
,Filling->Axis
,PlotStyle->{Thick,PointSize[Large]}
];
Show[Plot1,Plot2]

Section 1.4: Error

1.4.1: Definitions:

error[true_, approximate_]:=approximate-true
remainder[true_, approximate_]:=true - approximate
relativeError[true_, approximate_]:=error[true, approximate]/true

decimalsOfAccuracy[error_]:=-Log[10,Abs[error]]
digitsOfAccuracy[relativeError_]:=-Log[10,Abs[relativeError]]

eApprox=2.71828182 (* truncate e to eight decimals *)
eTrue=Exp[1]
N[eApprox,9]
N[eTrue,20]
error[Exp[1],eApprox]
relativeError[Exp[1],eApprox]
decimalsOfAccuracy[error[Exp[1],eApprox]]
digitsOfAccuracy[relativeError[Exp[1],eApprox]]

eApprox=2.71828183(* round e to eight decimals *)
error[Exp[1],eApprox]
decimalsOfAccuracy[error[Exp[1],eApprox]]
digitsOfAccuracy[relativeError[Exp[1],eTrue]]

1.4.2: Sources of error

f[x_]:=Exp[x]
fhat[x_]:=Sum[x^i/Factorial[i],{i,0,3}]
Plot[{f[x],fhat[x]},{x,0,2}]
a=Sqrt[2]
ahat=1.41
error[f[ahat],f[a]]
error[fhat[ahat],f[ahat]]