(* Illustration of Mathematica for simulating arbitrary precision arithmetic.*)

(* <<NumericalMath`ComputerArithmetic`  older Mathematica versions.*)
<<ComputerArithmetic`
SetArithmetic[5,10,MixedMode->True]  (* Sets 5-digit base 10 arithmetic.*)
x = {100.01, 102.01, .01}; y = {1.02, -1.00, .04};

(* Now convert each entry in these vectors to a computer number.*)
x5 = Map[ComputerNumber,x]; y5 = Map[ComputerNumber,y];

(* Here's the exact result, using the Mathematica default 64-bit arithmetic.*)
x . y ;

(* Here's the result, using 5-digit decimal arithmetic, implied by the
use of ComputerNumbers.*)
x5 . y5 ;  (* The semicolon inhibits echoing the result.*)
Clear[x,y,x5,y5]


(*Recursion.*)

Nest[f,t,3]
NestList[f,t,3]
g[t_] := (t^2 + 2)/(2*t)
NestList[g,1.0,5] (*This computes Sqrt[2].*)

(* Hint for Homework # 1.  Suppose you needed to compute a recursion of
the form   p_{n+1} = f[p_n, n] ; this is a bit subtle, because the
recursion depends not only on the previous value, but also a changing
index.*)

FoldList[f, p2, Range[2,5]]
Last[FoldList[f, p2, Range[2,5]]]
FoldList[Plus,0,{a,b,c}]  (*Just another example of FoldList.*)

(*To print out a specific number of digits, use NumberForm.*)
{N[Pi], NumberForm[N[Pi], 12]} (* N[] produces a numerical value.*)


(* Examples of function definition.*)

pnorm[x_?VectorQ, p_Real] := (Sum[(Abs[x[[i]]])^p,
	{i,1,Length[x]}])^(1/p) /; p >= 1.0

p[x_,limit_] := NSum[1/(k(k+x)),{k,limit,1,-1}, NSumTerms->limit]


(* Inner product examples.*)

ss = {{1, 0, 1}, {1, 0, 3}}; f = {1, 2, 3};
{\[Phi][0],\[Phi][1]} = Orthogonalize[ss, Dot]
(* Coefficients of f with respect to this orthogonal basis.*)
%.f
(* Projection of f onto subspace with basis ss.*)
(\[Phi][0].f) \[Phi][0] + (\[Phi][1].f) \[Phi][1]

(* Same examples with a function vector space and inner product
   (Integrate[#1 #2, {x, -1, 1}])& .*)
fss = {1, x, x^2, x^3}; fvec[x_] := Sin[x];
{\[Phi][0],\[Phi][1],\[Phi][2],\[Phi][3]} =
  Orthogonalize[fss, (Integrate[#1 #2, {x, -1, 1}]) &]
(* Projection of fvec onto the second orthonormal basis vector.*)
Projection[fvec[x], \[Phi][1], (NIntegrate[#1 #2, {x, -1, 1}]) &]
(* Best approximation to fvec by a polynomial of degree <= 3, with
   respect to the inner product (Integrate[#1 #2, {x, -1, 1}])& .*)
Sum[Projection[fvec[x], \[Phi][i],
  (NIntegrate[#1 #2, {x, -1, 1}])& ], {i, 0, 3}]


(* Examples of polynomial interpolation.*)

data={{0,12},{-1,-7},{2,20},{3,21}};
p[x_] = InterpolatingPolynomial[data,x]  (*Note use of = rather than := .*)
dataplot=ListPlot[data,PlotStyle->{PointSize[.03]}]
SetOptions[Plot, PlotRange->All];
pplot=Plot[p[x], {x,-1.1,3.1}]
Show[dataplot,pplot, PlotRange->All]

spline = Interpolation[data, InterpolationOrder->3]
splot=Plot[spline[x], {x,-1.0,3.0}, PlotStyle->{Dashing[{.05,.02}]}]
Show[dataplot, splot, PlotRange->All]


