Numerical Analysis studies algorithms which use numerical computations. The primary focus is on understanding errors. There are several types of errors:
Physical Measurement Errors. For example, the mass of the solar system is measured as \((1.98855 \pm 0.00025) \times 10^{30}\) kilograms. When you do calculations involving the mass of the solar system, you would do the entire calculation using \(1.98855\times 10^{30}\) but you would also perform some calculations to determine how big the error on the calculation is.
Floating point numbers have a limited number of bits so that all computations with floating point numbers are rounded. This rounding happens after every step of the computation. When doing a computation with floating point numbers, we want not only the result but also an estimate of the amount of error introduced by the rounding.
Mathematical approximation errors: For example, consider evaluating the integral \(\int_0^1 e^{-x^2} dx\). There is no anti-derivative for \(e^{-x^2}\) so the integral cannot be evaluated explicitly. Instead we would compute this using a Taylor series \[e^{-x^2} \approx 1 - x^2 + \frac{x^4}{2!} - \frac{x^6}{3!} + \frac{x^8}{4!}\] so we can approximate the integral \[\int_0^1 e^{-x^2} dx \approx \int_0^1 \left( 1 - x^2 + \frac{x^4}{2!} - \frac{x^6}{3!} + \frac{x^8}{4!} \right) dx\] And the above integral on the right can be computed with anti-derivatives. The error comes from the number of terms which were left out of the Taylor series approximation. Here as well we want an estimate for the error and would then like to keep track of how this error grows and propagates throughout future computation involving this integral.
One strange and interesting phenomenon is that if the same computation is performed in two different ways, the amount of error differs. For example, consider computing the following function on a \(4\)-decimal digit calculator (every computation is rounded to \(4\) decimal digits).
\[f(x) = x \left( \sqrt{x+1} - \sqrt{x} \right)\]
Consider \(x = 100\). Then \(\sqrt{100} = 10.0000\) and \(\sqrt{101} = 10.0499\) (note the rounding). Therefore the calculator computes
\[\sqrt{101} - \sqrt{100} = 10.0499 - 10.0000 = 0.0499\]
When multiplied by 100 we get
\[100 * 0.0499 = 4.9900\]
Notice the true value is \(4.98756..\). But now consider the following form of the function
\[f(x) = \frac{x}{\sqrt{x+1} + \sqrt{x}}\]
This is equivalent to the first form. (To see this, multiply the original \(f\) by \(\sqrt{x+1} + \sqrt{x}\) in both the numerator and denominator.) But now, computing again on a \(4\)-decimal digit calculator,
\[\sqrt{101} + \sqrt{100} = 20.0499\]
and
\[\frac{100}{20.0499} = 4.9875\]
much closer to the actual value of \(4.987562...\).
What happened is an instance of the following principle, called the loss of significance principle. When two nearly equal values are subtracted, leading significant digits are lost. We can see this in action because \(\sqrt{101} - \sqrt{100}\) had only three significant digits and this loss of significant digits impacted our future computations. The other computation did not have a corresponding loss of significant digits and so obtained a better result.
These types of issues come up a lot, especially when we use Taylor series (which is a lot).
\[e^{-5} = 1 + \frac{(-5)}{1!} + \frac{(-5)^2}{2!} + \frac{(-5)^3}{3!} + \frac{(-5)^4}{4!} \cdots\]
Note the terms are alternating between positive and negative so if we compute using this expansion, we will have lots of significance lost. If instead we use
\[e^{-5} = \frac{1}{e^5} = \frac{1}{\text{ series for } e^5}\]
the series for \(e^5\) has only positive terms so the approximation will be much better.
Thus changing the computation has an effect on the amount of accumulated error from roundoff error. The other types of error (measurement error and approximation errors) have the same feature so that computing the same thing in different ways impacts the amount of error that is propagated.
Numerical analysis in SAGE comes in several ways:
Some numerical methods attached to the symbolic classes we looked at last time. This allows you to mix symbolic and numeric computation. A big one is find_root. Example: FindRoot.sagews
Two python modules: numpy and scipy. These two modules are very popular for all areas of science. Numpy provides a very fast multidimensional array and Scipy contains a whole host of already implemented numerical algorithms. Since these are python modules, they can of course be used with SAGE.
A whole host of SAGE modules are built on top of numpy and scipy and provide even more advanced numerical analysis and of course the ability to combine numerical analysis with other Sage modules. (see here).
An example of the latter is towards the bottom of this article. Here is my sage worksheet: CurveFitting.sagews
find_root twice, once to find each intersection point since they intersect more than once. Why does solve not work?