A problem about fixed point iteration theory

Question

Description

Recently, I have been learning a couse called "Numerical Analysis". The fixed point iteration theory was introducted to solve the the approximate root of the equations in this course.

Fixed point iteration theory: $x=Ψ(x) \Rightarrow x_{k+1}=Ψ(x_k)$

Example

Solving the approximate root of $f(x)=x^3+4x^2-10 = 0$

So I give three styles of $x=Ψ(x)$, shown as below:

• $Ψ_1(x)=\frac{1}{2}\sqrt{10-x^3}$

• $Ψ_2(x)=\sqrt{\frac{10}{x+4}}$

• $Ψ_3(x)=x-\frac{f(x)}{f'(x)}=x-\frac{x^3+4x^2-10}{3x^2+8x}$

Using Mathematica to implement it and verify the results of textbook

My trial: I set the precision to 10

 result1 = FixedPointList[N[1/2 Sqrt[10 - #^3] &, 10], 1.5`10]

{1.500000000, 1.286953768, 1.402540804, 1.345458374, 1.375170253, 
   1.360094193, 1.367846968, 1.363887004, 1.365916733, 1.364878217, 
   1.365410061, 1.365137821, 1.365277209, 1.365205850, 1.365242384, 
   1.365223680, 1.365233256, 1.365228353, 1.365230863, 1.365229578, 
   1.365230236, 1.365229899, 1.365230072, 1.365229984, 1.365230029, 
   1.365230006, 1.365230017, 1.365230011, 1.365230014, 1.365230013, 
   1.365230014}

  result2 = FixedPointList[N[Sqrt[10/(# + 4)] &, 10], 1.5`10]

{1.500000000, 1.348399725, 1.367376372, 1.364957015, 1.365264748, 
   1.365225594, 1.365230576, 1.365229942, 1.365230023, 1.365230012, 
   1.365230014, 1.365230013}

 result3 = FixedPointList[N[# - (#^3 + 4 #^2 - 10)/(3 #^2 + 8 #) &, 10], 1.5`10]

 {1.500000000, 1.373333333, 1.365262015, 1.36523001, 1.36523001}

However, the last two results of the result3 is 1.36523001, 1.36523001, whose precision is 9 rather than 10

TableForm[
  Flatten[{result1, result2, result3}, {{2}, {1}}], 
  TableHeadings -> {None, Style[#, 15, Red] & /@ {"\[Phi]1 steps", "\[Phi]2 steps", 
   "\[Phi]3 steps"}}]

As the picture shown, the result2 ia same as the result of textbook, however, the result1 and result3 are a little different from the result of textbook. In addition, the precison of result3 is 9, not 10.

Question

Can someone give me a explalation about this trial? For me, I cannot understand this result.

Answer 1

A nice trick to force Mathematica to use a given precision is to use Block and make $MinPrecision equal to $MaxPrecision. So you can write your result1 as:

Block[{$MinPrecision = 10, $MaxPrecision = 10}, 
         FixedPointList[N[1/2 Sqrt[10 - #^3] &, 10], 1.5`10]]

{1.500000000, 1.286953768, 1.402540804, 1.345458374, 1.375170253, 1.360094193,
 1.367846968, 1.363887004, 1.365916733, 1.364878217, 1.365410061, 1.365137821,
 1.365277209, 1.365205850, 1.365242384, 1.365223680, 1.365233256, 1.365228353,
 1.365230863, 1.365229578,  1.365230236, 1.365229899, 1.365230072, 1.365229984,
 1.365230029, 1.365230006, 1.365230017, 1.365230011, 1.365230014, 1.365230013,      
 1.365230014, 1.365230013, 1.365230013, 1.365230013}

And similarly for result3

Block[{$MinPrecision = 10, $MaxPrecision = 10}, 
     FixedPointList[N[# - (#^3 + 4 #^2 - 10)/(3 #^2 + 8 #) &, 10], 1.5`10]]

{1.500000000, 1.373333333, 1.365262015, 1.365230014, 1.365230013, 1.365230013}

Answer 2

Mathematical operations on numbers of a given precision cannot guarantee the output to maintain the precision of the input numbers. In general precision will decrease. The amount of decrease depends on the operations and some operations will cause precision to decrease significantly. If you want a specific precision on the final result then the working precision should be higher than the desired output precision. Frequently, the working precision is done at twice the desired output precision to allow for significant lose of precision during complex calculations. In your example just using machine precision numbers will perform well.

result1 =
  FixedPointList[N[1/2 Sqrt[10 - #^3] &], 1.5];
{Length[result1], result1[[-1]] // InputForm}

{53, InputForm[1.3652300134140969`]}

result2 =
  FixedPointList[N[Sqrt[10/(# + 4)] &], 1.5];
{Length[result2], result2[[-1]] // InputForm}

{19, InputForm[1.3652300134140969`]}

result3 =
  FixedPointList[N[# - (#^3 + 4 #^2 - 10)/(3 #^2 + 8 #) &], 1.5];
{Length[result3], result3[[-1]] // InputForm}

{6, InputForm[1.3652300134140969`]}

You can shorten the displayed numbers using NumberForm