Jump to content

Numerical stability: Difference between revisions

From Wikipedia, the free encyclopedia
Content deleted Content added
mNo edit summary
fix Babylonian method formula
 
(46 intermediate revisions by 37 users not shown)
Line 1: Line 1:
{{Short description|Ability of numerical algorithms to remain accurate under small changes of inputs}}
{{No footnotes|date=February 2012}}
{{More footnotes|date=February 2012}}
In the [[mathematics|mathematical]] subfield of [[numerical analysis]], '''numerical stability''' is a desirable property of [[numerical algorithm]]s. The precise definition of stability depends on the context, but it is derived from the accuracy of the algorithm.
In the [[mathematics|mathematical]] subfield of [[numerical analysis]], '''numerical stability''' is a generally desirable property of [[numerical algorithm]]s. The precise definition of stability depends on the context. One is [[numerical linear algebra]] and the other is algorithms for solving ordinary and partial differential equations by discrete approximation.


In numerical linear algebra, the principal concern is instabilities caused by proximity to singularities of various kinds, such as very small or nearly colliding [[eigenvalues]]. On the other hand, in numerical algorithms for differential equations the concern is the growth of round-off errors and/or small fluctuations in initial data which might cause a large deviation of final answer from the exact solution.{{Citation needed|date=October 2017}}
An [[opposite (semantics)|opposite]] phenomenon is '''instability'''. Typically, algorithms would approach the right solution in the limit, if there were no round-off or truncation errors, but depending on the specific computational method, errors can be magnified, instead of damped, causing the error to grow exponentially.


Sometimes a single calculation can be achieved in several ways, all of which are algebraically equivalent in terms of ideal real or complex numbers, but in practice when performed on digital computers yield different results. Some calculations might damp out approximation errors that occur; others might magnify such errors. Calculations that can be proven not to magnify approximation errors are called ''numerically stable''. One of the common tasks of numerical analysis is to try to select algorithms which are ''robust'' – that is to say, have good numerical stability among other desirable properties.
Some numerical algorithms may damp out the small fluctuations (errors) in the input data; others might magnify such errors. Calculations that can be proven not to magnify approximation errors are called ''numerically stable''. One of the common tasks of numerical analysis is to try to select algorithms which are ''robust'' – that is to say, do not produce a wildly different result for very small change in the input data.


An [[opposite (semantics)|opposite]] phenomenon is '''instability'''. Typically, an algorithm involves an approximative method, and in some cases one could prove that the algorithm would approach the right solution in some limit (when using actual real numbers, not floating point numbers). Even in this case, there is no guarantee that it would converge to the correct solution, because the floating-point round-off or truncation errors can be magnified, instead of damped, causing the deviation from the exact solution to grow exponentially.<ref>{{cite book | title=Numerical Algorithms with C | author1=Giesela Engeln-Müllges|author2= Frank Uhlig |others= M. Schon (Translator), F. Uhlig (Translator) | edition=1 |date=2 July 1996 | publisher=Springer | pages=10 | isbn=978-3-540-60530-0 | url=https://books.google.com/books?id=HurESoDQljcC&pg=PA10}}</ref>
==Example==
As an example of an unstable algorithm, consider the task of adding an array of 100 numbers. To simplify things, assume our computer only has two [[significant figures|digits of precision]] (for example, numbers can be represented as 2.3, 77, 100, 110, 120, etc., but not 12.3 or 177).

The naive way to do this would be the following:
<code>
'''function''' sumArray(array) '''is'''
'''let''' theSum = 0
'''for each''' element '''in''' array '''do'''
'''let''' theSum = theSum + element
'''end for'''
'''return''' theSum
'''end function'''</code>

That looks reasonable, but suppose the first element in the array was 1.0 and the other 99 elements were 0.01. In exact arithmetic, the answer would be 1.99. However, on our two-digit computer, once the 1.0 was added into the sum variable, adding in 0.01 would have no effect on the sum, and so the final answer would be 1.0 – not a very good approximation of the real answer. Furthermore, we see that the algorithm depends on the ordering of elements within the array, in contrast to the exact arithmetic.


== Stability in numerical linear algebra ==
A stable algorithm would first sort the array by the absolute values of the elements in ascending order. This ensures that the numbers closest to zero will be taken into consideration first. Once that change is made, all of the 0.01 elements will be added, giving 0.99, and then the 1.0 element will be added, yielding a rounded result of 2.0 – a much better approximation of the real result.

==Forward, backward, and mixed stability==
There are different ways to formalize the concept of stability. The following definitions of forward, backward, and mixed stability are often used in [[numerical linear algebra]].
There are different ways to formalize the concept of stability. The following definitions of forward, backward, and mixed stability are often used in [[numerical linear algebra]].


[[Image:Forward and backward error.svg|frame|Diagram showing the forward error &Delta;''y'' and the backward error &Delta;''x'', and their relation to the exact solution map&nbsp;{{mvar|f}} and the numerical solution&nbsp;{{mvar|f}}*.]]
[[File:Forward and backward error.svg|frame|Diagram showing the '''forward error''' {{math|Δ''y''}} and the '''backward error''' {{math|Δ''x''}}, and their relation to the exact solution map&nbsp;{{mvar|f}} and the numerical solution&nbsp;{{mvar|f*}}.]]


Consider the problem to be solved by the numerical algorithm as a [[function (mathematics)|function]]&nbsp;{{mvar|f}} mapping the data&nbsp;{{mvar|x}} to the solution&nbsp;{{mvar|y}}. The result of the algorithm, say {{mvar|y}}*, will usually deviate from the "true" solution&nbsp;{{mvar|y}}. The main causes of error are [[round-off error]] and [[truncation error]]. The ''forward error'' of the algorithm is the difference between the result and the solution; in this case, {{math|1=Δ''y'' = ''y''* − ''y''}}. The ''backward error'' is the smallest Δ{{mvar|x}} such that {{math|1=''f'' (''x'' + Δ''x'') = ''y''*}}; in other words, the backward error tells us what problem the algorithm actually solved. The forward and backward error are related by the [[condition number]]: the forward error is at most as big in magnitude as the condition number multiplied by the magnitude of the backward error.
Consider the problem to be solved by the numerical algorithm as a [[function (mathematics)|function]]&nbsp;{{mvar|f}} mapping the data&nbsp;{{mvar|x}} to the solution&nbsp;{{mvar|y}}. The result of the algorithm, say {{mvar|y}}*, will usually deviate from the "true" solution&nbsp;{{mvar|y}}. The main causes of error are [[round-off error]] and [[truncation error]]. The ''forward error'' of the algorithm is the difference between the result and the solution; in this case, {{math|1=Δ''y'' = ''y''* − ''y''}}. The ''backward error'' is the smallest {{math|Δ''x''}} such that {{math|1=''f'' (''x'' + Δ''x'') = ''y''*}}; in other words, the backward error tells us what problem the algorithm actually solved. The forward and backward error are related by the [[condition number]]: the forward error is at most as big in magnitude as the condition number multiplied by the magnitude of the backward error.


In many cases, it is more natural to consider the [[relative error]]
In many cases, it is more natural to consider the [[relative error]] <math display="block"> \frac{|\Delta x|}{|x|} </math> instead of the absolute error {{math|Δ''x''}}.
:<math> \frac{|\Delta x|}{|x|} </math>
instead of the absolute error Δ{{mvar|x}}.


The algorithm is said to be ''backward stable'' if the backward error is small for all inputs&nbsp;{{mvar|x}}. Of course, "small" is a relative term and its definition will depend on the context. Often, we want the error to be of the same order as, or perhaps only a few [[orders of magnitude]] bigger than, the [[unit round-off]].
The algorithm is said to be ''backward stable'' if the backward error is small for all inputs&nbsp;{{mvar|x}}. Of course, "small" is a relative term and its definition will depend on the context. Often, we want the error to be of the same order as, or perhaps only a few [[orders of magnitude]] bigger than, the [[unit round-off]].


[[Image:Mixed stability diagram.svg|thumb|Mixed stability combines the concepts of forward error and backward error.]]
[[File:Mixed stability diagram.svg|thumb|Mixed stability combines the concepts of forward error and backward error.]]


The usual definition of numerical stability uses a more general concept, called ''mixed stability'', which combines the forward error and the backward error. An algorithm is stable in this sense if it solves a nearby problem approximately, i.e., if there exists a Δ{{mvar|x}} such that both Δ{{mvar|x}} is small and {{math|''f'' (''x'' + Δ''x'') − ''y''*}} is small. Hence, a backward stable algorithm is always stable.
The usual definition of numerical stability uses a more general concept, called ''mixed stability'', which combines the forward error and the backward error. An algorithm is stable in this sense if it solves a nearby problem approximately, i.e., if there exists a {{math|Δ''x''}} such that both {{math|Δ''x''}} is small and {{math|''f'' (''x'' + Δ''x'') − ''y''*}} is small. Hence, a backward stable algorithm is always stable.


An algorithm is ''forward stable'' if its forward error divided by the condition number of the problem is small. This means that an algorithm is forward stable if it has a forward error of magnitude similar to some backward stable algorithm.
An algorithm is ''forward stable'' if its forward error divided by the condition number of the problem is small. This means that an algorithm is forward stable if it has a forward error of magnitude similar to some backward stable algorithm.


== Stability in numerical differential equations ==
==Error growth==
The above definitions are particularly relevant in situations where truncation errors are not important. In other contexts, for instance when solving [[differential equation]]s, a different definition of numerical stability is used.
[[Image:errorgrowth.jpg|300px|thumb|right|Comparing the linear error growth of a stable algorithm and the exponential error growth of an unstable algorithm used to solve the same problem, with the same initial data.]]


In [[Numerical methods for ordinary differential equations|numerical ordinary differential equations]], various concepts of numerical stability exist, for instance [[Stiff equation#A-stability|A-stability]]. They are related to some concept of stability in the [[dynamical system]]s sense, often [[Lyapunov stability]]. It is important to use a stable method when solving a [[stiff equation]].
Suppose that {{math|''E''<sub>''i''</sub> > 0}} denotes an initial error and {{math|''E''<sub>''n''</sub>}} represents the magnitude of an error after {{mvar|n}} subsequent operations. If {{math|''E''<sub>''n''</sub> ∼  ''C''∙''n''∙''E''<sub>''i''</sub>}}, where {{mvar|C}} is a constant independent of&nbsp;{{mvar|n}}, then the [[asymptotic analysis|growth]] of the error is said to be ''linear''. If {{math|''E''<sub>''n''</sub> ∼  ''C''<sup>''n''</sup>∙''E''<sub>''i''</sub>}}, for some {{math|''C'' > 1}}, then the growth of the error is called [[exponential growth|''exponential'']].


Yet another definition is used in [[numerical partial differential equations]]. An algorithm for solving a linear evolutionary [[partial differential equation]] is stable if the [[total variation]] of the numerical solution at a fixed time remains bounded as the step size goes to zero. The [[Lax equivalence theorem]] states that an algorithm [[Numerical methods for ordinary differential equations#Convergence|converges]] if it is [[Numerical methods for ordinary differential equations#Consistency and order|consistent]] and [[Numerical methods for ordinary differential equations#Stability and stiffness|stable]] (in this sense). Stability is sometimes achieved by including [[numerical diffusion]]. Numerical diffusion is a mathematical term which ensures that roundoff and other errors in the calculation get spread out and do not add up to cause the calculation to "blow up". [[Von Neumann stability analysis]] is a commonly used procedure for the stability analysis of [[finite difference method|finite difference schemes]] as applied to linear partial differential equations. These results do not hold for nonlinear PDEs, where a general, consistent definition of stability is complicated by many properties absent in linear equations.
==Stability in numerical differential equations==
The above definitions are particularly relevant in situations where truncation errors are not important. In other contexts, for instance when solving [[differential equation]]s, a different definition of numerical stability is used.


==Example==
In [[numerical ordinary differential equations]], various concepts of numerical stability exist, for instance [[A-stability]]. They are related to some concept of stability in the [[dynamical systems]] sense, often [[Lyapunov stability]]. It is important to use a stable method when solving a [[stiff equation]].
Computing the square root of 2 (which is roughly 1.41421) is a [[well-posed problem]]. Many algorithms solve this problem by starting with an initial approximation ''x''<sub>0</sub> to <math>\sqrt{2}</math>, for instance ''x''<sub>0</sub> = 1.4, and then computing improved guesses ''x''<sub>1</sub>, ''x''<sub>2</sub>, etc. One such method is the famous [[Babylonian method]], which is given by ''x''<sub>''k''+1</sub> = (''x<sub>k</sub>''+ 2/''x<sub>k</sub>'')/2. Another method, called "method X", is given by ''x''<sub>''k''+1</sub> = (''x''<sub>''k''</sub><sup>2</sup> − 2)<sup>2</sup> + ''x''<sub>''k''</sub>.{{NoteTag|This is a [[fixed point iteration]] for the equation <math>x=(x^2-2)^2+x=f(x)</math>, whose solutions include <math>\sqrt{2}</math>. The iterates always move to the right since <math>f(x)\geq x</math>. Hence <math>x_1=1.4<\sqrt{2}</math> converges and <math>x_1=1.42>\sqrt{2}</math> diverges.}} A few iterations of each scheme are calculated in table form below, with initial guesses ''x''<sub>0</sub> = 1.4 and ''x''<sub>0</sub> = 1.42.


{| class="wikitable"
Yet another definition is used in [[numerical partial differential equations]]. An algorithm for solving a linear evolutionary [[partial differential equation]] is stable if the [[total variation]] of the numerical solution at a fixed time remains bounded as the step size goes to zero. The [[Lax equivalence theorem]] states that an algorithm converges if it is consistent and stable (in this sense). Stability is sometimes achieved by including [[numerical diffusion]]. Numerical diffusion is a mathematical term which ensures that roundoff and other errors in the calculation get spread out and do not add up to cause the calculation to "blow up". [[von Neumann stability analysis]] is a commonly used procedure for the stability analysis of [[finite difference scheme]]s as applied to linear partial differential equations. These results do not hold for nonlinear PDEs, where a general, consistent definition of stability is complicated by many properties absent in linear equations.
|-
! Babylonian
! Babylonian
! Method X
! Method X
|-
| ''x''<sub>0</sub> = 1.4
| ''x''<sub>0</sub> = 1.42
| ''x''<sub>0</sub> = 1.4
| ''x''<sub>0</sub> = 1.42
|-
| ''x''<sub>1</sub> = 1.4142857...
| ''x''<sub>1</sub> = 1.41422535...
| ''x''<sub>1</sub> = 1.4016
| ''x''<sub>1</sub> = 1.42026896
|-
| ''x''<sub>2</sub> = 1.414213564...
| ''x''<sub>2</sub> = 1.41421356242...
| ''x''<sub>2</sub> = 1.4028614...
| ''x''<sub>2</sub> = 1.42056...
|-
|
|
| ...
| ...
|-
|
|
| ''x''<sub>1000000</sub> = 1.41421...
| ''x''<sub>27</sub> = 7280.2284...
|}


Observe that the Babylonian method converges quickly regardless of the initial guess, whereas Method X converges extremely slowly with initial guess ''x''<sub>0</sub> = 1.4 and diverges for initial guess ''x''<sub>0</sub> = 1.42. Hence, the Babylonian method is numerically stable, while Method X is numerically unstable.
==See also==

Numerical stability is affected by the number of the significant digits the machine keeps. If a machine is used that keeps only the four most significant decimal digits, a good example on loss of significance can be given by the two equivalent functions
:<math>
f(x)=x\left(\sqrt{x+1}-\sqrt{x}\right)
</math> and <math>
g(x)=\frac{x}{\sqrt{x+1}+\sqrt{x}}.
</math>
:Comparing the results of
:: <math> f(500)=500 \left(\sqrt{501}-\sqrt{500} \right)=500 \left(22.38-22.36 \right)=500(0.02)=10</math>
:and
:<math>
\begin{alignat}{3}g(500)&=\frac{500}{\sqrt{501}+\sqrt{500}}\\
&=\frac{500}{22.38+22.36}\\
&=\frac{500}{44.74}=11.17
\end{alignat}
</math>
by comparing the two results above, it is clear that [[loss of significance]] (caused here by [[catastrophic cancellation]] from subtracting approximations to the nearby numbers <math>\sqrt{501}</math> and <math>\sqrt{500}</math>, despite the subtraction being computed exactly) has a huge effect on the results, even though both functions are equivalent, as shown below
:<math> \begin{alignat}{4}
f(x)&=x \left(\sqrt{x+1}-\sqrt{x} \right)\\
&=x \left(\sqrt{x+1}-\sqrt{x} \right)\frac{\sqrt{x+1}+\sqrt{x}}{\sqrt{x+1}+\sqrt{x}}\\
&=x\frac{(\sqrt{x+1})^2-(\sqrt{x})^2}{\sqrt{x+1}+\sqrt{x}}\\
&=x\frac{x+1-x}{\sqrt{x+1}+\sqrt{x}} \\
&=x\frac{1}{\sqrt{x+1}+\sqrt{x}} \\
&=\frac {x}{\sqrt{x+1}+\sqrt{x}} \\
&=g(x)
\end{alignat}</math>

The desired value, computed using infinite precision, is 11.174755...{{NoteTag|The example is a modification of one taken from {{harvtxt|Mathews|Fink|1999}}.<ref>{{cite book|page=28|contribution=Example 1.17|title=Numerical Methods Using MATLAB|edition=3rd|first1=John H.|last1=Mathews|first2=Kurtis D.|last2=Fink|publisher=Prentice Hall|year=1999}}</ref>}}

== See also ==
* [[Algorithms for calculating variance]]
* [[Algorithms for calculating variance]]
* [[Stability theory]]
* [[Chaos theory]]
* [[Chaos theory]]
* [[Propagation of uncertainty]]


==References==
==Notes==
{{reflist|group=note}}
*[[Nicholas J. Higham]], ''Accuracy and Stability of Numerical Algorithms'', Society of Industrial and Applied Mathematics, Philadelphia, 1996. ISBN 0-89871-355-2.
*Richard L. Burden and J. Douglas Faires, ''Numerical Analysis 8th Edition'', Thomson Brooks/Cole, U.S., 2005. ISBN 0-534-39200-8


== References ==
{{DEFAULTSORT:Numerical Stability}}
{{Reflist}}
[[Category:Numerical analysis]]
* {{cite book|author= Nicholas J. Higham|author-link= Nicholas Higham|title= Accuracy and Stability of Numerical Algorithms|publisher= Society of Industrial and Applied Mathematics|location= Philadelphia|date= 1996|isbn= 0-89871-355-2|url-access= registration|url= https://archive.org/details/accuracystabilit0000high}}
* {{cite book|author1=Richard L. Burden |author2= J. Douglas Faires|title=Numerical Analysis|edition= 8th|publisher= Thomson Brooks/Cole|location= U.S.|date= 2005|isbn= 0-534-39200-8}}
* {{cite journal|arxiv=1605.04339 | doi=10.1109/MCSE.2017.3151254 |last1=Mesnard |first1=Olivier |title= Reproducible and Replicable Computational Fluid Dynamics: It's Harder Than You Think|last2= Barba |first2=Lorena A. | journal=Computing in Science & Engineering | year=2017 | volume=19 | issue=4 | pages=44–55 | bibcode=2017CSE....19d..44M | s2cid=11288122 }}


{{Linear algebra}}
[[ar:استقرار عددي]]

[[cs:Stabilita numerické metody]]
[[Category:Numerical analysis]]
[[de:Stabilität (Numerik)]]
[[es:Estabilidad numérica]]
[[eo:Cifereca stabileco]]
[[fr:Stabilité numérique]]
[[it:Stabilità numerica]]
[[ja:数値的安定性]]
[[pl:Algorytm numerycznie stabilny]]
[[pt:Estabilidade numérica]]
[[ru:Вычислительная устойчивость]]
[[ur:عددی ثبات]]
[[zh:数值稳定性]]

Latest revision as of 02:37, 26 February 2024

In the mathematical subfield of numerical analysis, numerical stability is a generally desirable property of numerical algorithms. The precise definition of stability depends on the context. One is numerical linear algebra and the other is algorithms for solving ordinary and partial differential equations by discrete approximation.

In numerical linear algebra, the principal concern is instabilities caused by proximity to singularities of various kinds, such as very small or nearly colliding eigenvalues. On the other hand, in numerical algorithms for differential equations the concern is the growth of round-off errors and/or small fluctuations in initial data which might cause a large deviation of final answer from the exact solution.[citation needed]

Some numerical algorithms may damp out the small fluctuations (errors) in the input data; others might magnify such errors. Calculations that can be proven not to magnify approximation errors are called numerically stable. One of the common tasks of numerical analysis is to try to select algorithms which are robust – that is to say, do not produce a wildly different result for very small change in the input data.

An opposite phenomenon is instability. Typically, an algorithm involves an approximative method, and in some cases one could prove that the algorithm would approach the right solution in some limit (when using actual real numbers, not floating point numbers). Even in this case, there is no guarantee that it would converge to the correct solution, because the floating-point round-off or truncation errors can be magnified, instead of damped, causing the deviation from the exact solution to grow exponentially.[1]

Stability in numerical linear algebra

[edit]

There are different ways to formalize the concept of stability. The following definitions of forward, backward, and mixed stability are often used in numerical linear algebra.

Diagram showing the forward error Δy and the backward error Δx, and their relation to the exact solution map f and the numerical solution f*.

Consider the problem to be solved by the numerical algorithm as a function f mapping the data x to the solution y. The result of the algorithm, say y*, will usually deviate from the "true" solution y. The main causes of error are round-off error and truncation error. The forward error of the algorithm is the difference between the result and the solution; in this case, Δy = y* − y. The backward error is the smallest Δx such that f (x + Δx) = y*; in other words, the backward error tells us what problem the algorithm actually solved. The forward and backward error are related by the condition number: the forward error is at most as big in magnitude as the condition number multiplied by the magnitude of the backward error.

In many cases, it is more natural to consider the relative error instead of the absolute error Δx.

The algorithm is said to be backward stable if the backward error is small for all inputs x. Of course, "small" is a relative term and its definition will depend on the context. Often, we want the error to be of the same order as, or perhaps only a few orders of magnitude bigger than, the unit round-off.

Mixed stability combines the concepts of forward error and backward error.

The usual definition of numerical stability uses a more general concept, called mixed stability, which combines the forward error and the backward error. An algorithm is stable in this sense if it solves a nearby problem approximately, i.e., if there exists a Δx such that both Δx is small and f (x + Δx) − y* is small. Hence, a backward stable algorithm is always stable.

An algorithm is forward stable if its forward error divided by the condition number of the problem is small. This means that an algorithm is forward stable if it has a forward error of magnitude similar to some backward stable algorithm.

Stability in numerical differential equations

[edit]

The above definitions are particularly relevant in situations where truncation errors are not important. In other contexts, for instance when solving differential equations, a different definition of numerical stability is used.

In numerical ordinary differential equations, various concepts of numerical stability exist, for instance A-stability. They are related to some concept of stability in the dynamical systems sense, often Lyapunov stability. It is important to use a stable method when solving a stiff equation.

Yet another definition is used in numerical partial differential equations. An algorithm for solving a linear evolutionary partial differential equation is stable if the total variation of the numerical solution at a fixed time remains bounded as the step size goes to zero. The Lax equivalence theorem states that an algorithm converges if it is consistent and stable (in this sense). Stability is sometimes achieved by including numerical diffusion. Numerical diffusion is a mathematical term which ensures that roundoff and other errors in the calculation get spread out and do not add up to cause the calculation to "blow up". Von Neumann stability analysis is a commonly used procedure for the stability analysis of finite difference schemes as applied to linear partial differential equations. These results do not hold for nonlinear PDEs, where a general, consistent definition of stability is complicated by many properties absent in linear equations.

Example

[edit]

Computing the square root of 2 (which is roughly 1.41421) is a well-posed problem. Many algorithms solve this problem by starting with an initial approximation x0 to , for instance x0 = 1.4, and then computing improved guesses x1, x2, etc. One such method is the famous Babylonian method, which is given by xk+1 = (xk+ 2/xk)/2. Another method, called "method X", is given by xk+1 = (xk2 − 2)2 + xk.[note 1] A few iterations of each scheme are calculated in table form below, with initial guesses x0 = 1.4 and x0 = 1.42.

Babylonian Babylonian Method X Method X
x0 = 1.4 x0 = 1.42 x0 = 1.4 x0 = 1.42
x1 = 1.4142857... x1 = 1.41422535... x1 = 1.4016 x1 = 1.42026896
x2 = 1.414213564... x2 = 1.41421356242... x2 = 1.4028614... x2 = 1.42056...
... ...
x1000000 = 1.41421... x27 = 7280.2284...

Observe that the Babylonian method converges quickly regardless of the initial guess, whereas Method X converges extremely slowly with initial guess x0 = 1.4 and diverges for initial guess x0 = 1.42. Hence, the Babylonian method is numerically stable, while Method X is numerically unstable.

Numerical stability is affected by the number of the significant digits the machine keeps. If a machine is used that keeps only the four most significant decimal digits, a good example on loss of significance can be given by the two equivalent functions

and
Comparing the results of
and

by comparing the two results above, it is clear that loss of significance (caused here by catastrophic cancellation from subtracting approximations to the nearby numbers and , despite the subtraction being computed exactly) has a huge effect on the results, even though both functions are equivalent, as shown below

The desired value, computed using infinite precision, is 11.174755...[note 2]

See also

[edit]

Notes

[edit]
  1. ^ This is a fixed point iteration for the equation , whose solutions include . The iterates always move to the right since . Hence converges and diverges.
  2. ^ The example is a modification of one taken from Mathews & Fink (1999).[2]

References

[edit]
  1. ^ Giesela Engeln-Müllges; Frank Uhlig (2 July 1996). Numerical Algorithms with C. M. Schon (Translator), F. Uhlig (Translator) (1 ed.). Springer. p. 10. ISBN 978-3-540-60530-0.
  2. ^ Mathews, John H.; Fink, Kurtis D. (1999). "Example 1.17". Numerical Methods Using MATLAB (3rd ed.). Prentice Hall. p. 28.