Archive for the ‘Numerical Methods’ Category

Newton’s method using rational arithmetic

April 2, 2012

The floating-point number system puts further restraints on what can be expected. We can at best hope to keep the relative error under control, and we cannot expect to find zeros far from the origin with great absolute accuracy.

– Richard W. Hamming [1]

I am somewhat disenchanted with floating-point arithmetic. We all have implemented numerical methods using floating-point numbers. But, how many of us have implemented numerical methods using rational numbers of arbitrary precision? I certainly have not. Hence, let us give it a try. In this post I will implement Newton’s method [1] (also known as the Newton-Raphson method) in Haskell using arbitrary-precision rational numbers.

Given a continuous (real) function f : \mathbb{R} \to \mathbb{R} and x_0, an initial approximation of the zero of function f, we then obtain the next approximation using the following recurrence relation

x_{k+1} = x_k - \displaystyle\frac{f (x_k)}{f' (x_k)}

where f' is the first derivative of f. We will use arbitrary-precision rational numbers (represented as a ratio of two Integer values)

type Rational = Ratio Integer

which are defined in the Data.Ratio library.

__________

Example

Suppose we would like to find a rational approximation of \sqrt{2}. We create a real function defined by f (x) = x^2 - 2, whose zeros are  \sqrt{2} and -\sqrt{2}. Note that the first derivative of f is f' (x) = 2 x. The initial approximation is, say, x_0 = 1 / 1, which is closer to \sqrt{2} than to -\sqrt{2}. Hopefully, the sequence \{x_k\} will converge to \sqrt{2} instead of converging to -\sqrt{2}. We are performing a numerical experiment, and success is not guaranteed!

The following Haskell script implements Newton’s method:

import Data.Ratio

-- define function f
f :: Rational -> Rational
f x = x*x - 2

-- define 1st derivative of function f
f' :: Rational -> Rational
f' x = 2*x

-- define Newton map
g :: Rational -> Rational
g x = x - (f x / f' x)

-- initial approximation
x0 :: Rational
x0 = 1 % 1

-- list of rational approximations
xs :: [Rational]
xs = iterate g x0

where we used (higher-order) function iterate to create the list of rational approximations. We run the script above in GHCi and use take to compute the first few rational approximations in the list:

*Main> take 6 xs
[1 % 1,3 % 2,17 % 12,577 % 408,665857 % 470832,
886731088897 % 627013566048]

Hence, we have, for instance, the following rational approximation

x_5 = \displaystyle\frac{886731088897}{627013566048} \approx \sqrt{2}.

How good is this rational approximation? Let us see:

*Main> let x5 = xs !! 5
*Main> x5*x5 - 2
1 % 393146012008229658338304
*Main> 886731088897 / 627013566048
1.4142135623730951
*Main> sqrt 2
1.4142135623730951

Fairly good indeed! And all it took was five iterations!

__________

References

[1] Richard W. Hamming, Numerical Methods for Scientists and Engineers, Dover Publications, 1973.

Tags:, , , , , , ,
Posted in Haskell, Numerical Methods | 10 Comments »

Interpolation and energy amplification

March 25, 2011

Consider the following interpolator [1]

where x_u = (\uparrow L) (x) is the upsampling of sequence x by an integer factor of L, and where H(z) is the transfer function of an ideal discrete-time LTI lowpass filter with gain L and cutoff frequency \pi/L. What is the relation between the energies of signals x and y?

__________

My solution:

It can be easily shown [1] that the interpolator’s input-output relationship is Y (z) = H(z) X (z^L). Hence, evaluating the Z-transforms on the unit circle, we obtain

Y (e^{j \omega}) = H(e^{j \omega}) X (e^{j \omega L}).

The energy of the interpolator’s output signal is

\mathcal{E}_y := \displaystyle \sum_{n \in \mathbb{Z}} | y(n) |^2 = \frac{1}{2 \pi} \int_{-\pi}^{\pi} |Y (e^{j \omega})|^2 d \omega

where we used Parseval’s theorem for discrete-time signals. Thus,

\mathcal{E}_y = \displaystyle\frac{1}{2 \pi} \int_{-\pi}^{\pi} |H(e^{j \omega})|^2 |X (e^{j \omega L})|^2 d \omega.

H(z) is an ideal lowpass filter with the following frequency response

H(e^{j \omega}) = \begin{cases} L, & |\omega| \leq \pi/L\\ 0, & \pi/ L< |\omega| \leq \pi\end{cases}

and the following sinc impulse response

h(n) = \displaystyle\frac{1}{2 \pi} \int_{-\pi}^{\pi} H(e^{j \omega}) e^{j \omega n} d \omega = \displaystyle \frac{\sin(\frac{n \pi}{L})}{\frac{n \pi}{L}} =: \text{sinc} \left(\frac{n \pi}{L}\right).

Please do note that this is a non-causal filter with an impulse response of infinite length. We can then write the interpolator’s output signal’s energy as follows

\mathcal{E}_y = \displaystyle\frac{1}{2 \pi} \int_{-\pi}^{\pi} |H(e^{j \omega})|^2 |X (e^{j \omega L})|^2 d \omega = \displaystyle\frac{L^2}{2 \pi} \int_{-\pi/L}^{\pi/L} |X (e^{j \omega L})|^2 d \omega.

Performing a change of variable, \theta = L \omega, we obtain

\mathcal{E}_y = \displaystyle\frac{L^2}{2 \pi} \int_{-\pi/L}^{\pi/L} |X (e^{j \omega L})|^2 d \omega = \displaystyle\frac{L}{2 \pi} \int_{-\pi}^{\pi} |X (e^{j \theta})|^2 d \theta = L \mathcal{E}_x

where

\mathcal{E}_x := \displaystyle \sum_{n \in \mathbb{Z}} | x(n) |^2 = \frac{1}{2 \pi} \int_{-\pi}^{\pi} |X (e^{j \omega})|^2 d \omega

is the energy of the input signal. Since \mathcal{E}_y = L \mathcal{E}_x, we can conclude that the interpolator amplifies the input signal’s energy by a factor of L, which is somewhat intuitive.

We considered the case where signals are bi-infinite sequences and H(z) is a non-causal ideal lowpass filter. It would be interesting to consider the more realistic case where the signals are finite sequences and H(z) is a causal FIR filter.

__________

References

[1] Alan V. Oppenheim, Ronald W. Schafer, Discrete-Time Signal Processing, 2nd edition, Prentice-Hall, 1999.

Tags:, , , ,
Posted in Numerical Methods, Puzzles, Signal Processing | 2 Comments »

Variational Integrators

June 3, 2008

In my early days as an undergraduate student, I had a lot of fun writing C code to integrate ODE’s numerically, so I could simulate physical systems in a computer (integrating the pendulum ODE is a classic). However, I was not very happy to realize that traditional numerical methods very easily give results with no physical meaning.

For example, think of a simple pendulum oscillating on a vertical plane. Assuming that there are no damping forces, there’s energy conservation. However, traditional numerical methods (such as Runge-Kutta methods) do not, in general, preserve energy. The numerical simulation of the pendulum may thus exhibit non-physical energy dissipation (or amplification).

[ image courtesy of Ari Stern ]

When we integrate an ODE numerically, we obtain a discrete trajectory which we hope is “close enough” to the exact flow of the differential equation. In some problems in Physics, the continuous flow lies on a manifold, and we would like to preserve this geometric property: it would be wonderful if the discrete trajectory lied on the same manifold as the continuous flow. While thinking about numerical methods that preserve certain invariants, I became interested in the field of geometric integration.

One particular class of geometric integrators I have been interested in over the past year is the class of variational integrators. Many physical laws are formulated as variational principles from which differential equations can be derived. Instead of discretizing the differential equations, why not discretize the variational principles instead? This rather sexy thought lies at the basis of the nascent field of variational integrators.

Here’s an overview on variational integrators, by Prof. Matthew West. I have been reading some papers on the topic. This field seems to have been explored rather vigorously at Caltech, as there are/were a lot of people working on it:

and their collaborators. I have compiled a list of lecture slides, theses and papers on the area of variational integrators. Hope you find the list useful (I sure do!).

-/-

Some lecture slides for a gentle introduction:

-/-

A few PhD theses on variational integrators and their application:

-/-

Some papers on the foundations of variational integrators:

Applications to Computational Mechanics:

Applications to Optimal Control Theory:

Applications to Computational Electromagnetism:

Applications to Computational Physics:

The papers I have listed so far refer to deterministic integrators. However, quite recently Prof. Houman Owhadi and Nawaf Bou-Rabee have been working on stochastic variational integrators:

I wonder if stochastic variational integrators could be used in quantitative finance. Just a wild thought…

-/-

If you know some more papers, please let me know! Last but not least, Ganesh Swami once wrote a brief post on variational integrators.

Tags:, , , , , , , , , , , , , , , ,
Posted in Mathematics, Numerical Methods | 1 Comment »

Follow

Get every new post delivered to your Inbox.

Join 47 other followers