Math

DNA and Bessel functions

Posted on by John

I was reading a book on the history of the discovery of the structure of DNA [1] and was surprised by a few passing references to Bessel functions.

According to Claude,

When X-rays are diffracted by a helical structure, the resulting diffraction pattern breaks into a series of horizontal “layer lines.” Cochran, Crick, and Vand showed mathematically that the diffracted amplitude on the n-th layer line is proportional to a Bessel function of the first kind, order n:

Jn(2π r R)

where r is the radius of the helix and R is the distance out from the center (the meridian) in the diffraction pattern.

The citation for this paragraph is a paper from 1952 [2] that amazingly is behind a paywall.

Update: Evan pointed out something fishy in the expression above: “Isn’t there a dimensional problem with J(2π r R)? Maybe 2πr/R instead?”

This was Claude’s response:

  • r is the radius of the helix, a real-space length (say, in Å).
  • R is the radial coordinate in reciprocal space — it has units of inverse length (Å⁻¹), because that’s how X-ray diffraction/Fourier space works (spatial frequency, not distance).

So Evan and Claude were both right. There is an implicit division happening, but crystallographers do write r R.

Related posts

[1] Watson and Crick didn’t “discover DNA” as is commonly said. DNA was discovered in 1878. Watson and Crick discovered the structure of DNA in 1953.

[2] Cochran, W., Crick, F. H. C., & Vand, V. (1952). “The Structure of Synthetic Polypeptides. I. The Transform of Atoms on a Helix.” Acta Crystallographica, 5(5), 581–586.

Aitken acceleration before Aitken

Posted on by John

Kepler solved his eponymous equation

M = E − e sin(E)

by finding a fixed point of

E = M + e sin(E).

So guess a value of E and stick it into the right hand side. Then plug that value into the right hand side again. Kepler said a couple iterations should be enough. And a couple iterations are enough if the eccentricity e is small and you don’t need much accuracy.

The rate of convergence is determined by e. Kepler implicitly had in mind small values of e because he wasn’t aware of anything orbiting the sun in a highly elliptical orbit. Here’s an example with eccentricity 0.05, about the eccentricity of the orbits of Jupiter and Saturn.

from math import sin

M, E, e = 1, 1, 0.05
for _ in range(5):
    E = M + e*sin(E)
    residual = M - (E - e*sin(E))
    print(residual)

The residual after just two iterations is 2.77 × 10−5. If you change e to 0.2, the eccentricity of Mercury’s orbit, it takes three iterations to get comparable accuracy. Mercury has the most eccentric orbit of any object Kepler would have known about.

Now suppose you’d like to solve for E when M = 1 for Halley’s comet, and you’d like an error of less than 10−8. Now you need 16 iterations.

C. F. W. Peters discovered a faster algorithm in 1891.

E1 = M + e sin(E0)
E2 = M + e sin(E1)
E3 = (E2 E0 − E1²)/(E2 − 2E1 + E0)

Let’s look at the results of doing three iterations of Peters’ method for Halley’s comet.

M, E0, e = 1, 1, 0.967
for _ in range(3):
    E1 = M + e*sin(E0)
    E2 = M + e*sin(E1)
    E3 = (E2*E0 - E1**2)/(E2 - 2*E1 + E0)
    residual = M - (E - e*sin(E3))
    print(residual)
    E0 = E3 # for next iteration

This gives a residual of −7.23 × 10−10. Each iteration of Peters’ method requires a little more than twice as much work as an iteration of Kepler’s method, but 3 iterations of Peters’ method accomplished more than 16 iterations of Kepler’s method.

Peters’ algorithm from 1891 was a special case of Alexander Aitken’s series acceleration method published in 1926.

Related posts

Tighter bounds on alternating series remainder

Posted on by John

The alternating series test is part of the standard calculus curriculum. It says that if you truncate an alternating series, the remainder is bounded by the first term that was left out. This fact goes by in a blur for most students, but it becomes useful later if you need to do numerical computing.

To be more precise, assume we have a series of the form

\sum_{i=1}^\infty (-1)^i a_i

where the ai are positive and monotonically converge to zero. Then the tail of the series is bounded by its first term:

\left|R_n\right| = \left| \sum_{i=n+1}^\infty (-1)^i a_i \right| \leq a_{n+1}

The more we can say about the behavior of the ai the more we can say about the remainder. So far we’ve assumed that these terms go monotonically to zero. If their differences

\Delta a_i = a_i - a_{i+1}

also go monotonically to zero, then we have an upper and lower bound on the truncation error:

\frac{a_{n+1}}{2} \leq |R_n| \leq \frac{a_n}{2}

If the differences of the differences,

\Delta^2 a_i = \Delta (\Delta a_i)

also converge monotonically to zero, we can get a larger lower bound and a smaller upper bound on the remainder. In general, if the differences up to order k of the ai go to zero monotonically, then the remainder term can be bounded as follows.

\frac{a_{n+1}}{2} +\frac{\Delta a_{n+1}}{2^2} +\cdots+ \frac{\Delta^k a_{n+1}}{2^{k+1}} < \left|R_n\right| < \frac{a_n}{2} -\left\{ \frac{\Delta a_n}{2^2} +\cdots+ \frac{\Delta^k a_n}{2^{k+1}} \right\}.

Source: Mark B. Villarino. The Error in an Alternating Series. American Mathematical Monthly, April 2018, pp. 360–364.

Related posts

Golden iteration

Posted on by John

The expression

\sqrt{1 + \sqrt{1 + \sqrt{1 + \cdots}}}

converges to the golden ratio φ. Another way to say this is that the sequence defined by x0 = 1 and

x_{n+1} = \sqrt{x_n + 1}

for n > 0 converges to φ. This post will be about how it converges.

I wrote a little script to look at the error in approximating φ by xn and noticed that the error is about three times smaller at each step. Here’s why that observation was correct.

The ratio of the error at one step to the error at the previous step is

\frac{\sqrt{x + 1} - \varphi}{x - \varphi}

If x = φ + ε the expression above becomes

\frac{\sqrt{\varphi + \epsilon + 1} - \varphi}{\epsilon} = \frac{1}{\sqrt{2(3 + \sqrt{5})}} - \frac{\epsilon}{\sqrt{8} (3 + \sqrt{5})^{3/2}} + {\cal O}(\epsilon^2)

when you expand as a Taylor series in ε centered at 0. This says the error multiplied by a factor of about

\frac{1}{\sqrt{2(3 + \sqrt{5})}} = 0.309016994\ldots

at each step. The next term in the Taylor series is approximately −0.03ε, so the exact rate of convergence is a slightly faster at first, but essentially the error is multiplied by 0.309 at each iteration.

Solving H_n = 100

Posted on by John

The previous post includes code for solving the equation

Hn = m

i.e. finding the value of n for which the nth harmonic number is the closest to m. It works well for small values of m. It works for large m in the sense that the solution is very close to m, but it’s not necessarily the best solution.

For example, set m = 100. The code returns

n = 15092688622113830917200248731913020965388288

and indeed for that value of n,

Hn − 100 ≈ 3 × 10−15

and that’s as much as we could hope for with IEEE 754 floats.

The approximation

n = exp(m −γ)

is very good for large values of m. Using Mathematica we can find the exact value of n.

f[n_] := Log[n] + EulerGamma + 1/(2 n) - 1/(12 n^2)
n = Floor[Exp[100 - EulerGamma]];
N[f[n], 50]
100.00000000000000000000000000000000000000000000900
N[f[n - 1], 50]
99.999999999999999999999999999999999999999999942747

So

n = 15092688622113788323693563264538101449859497

A similar process can find the solution to

Hn = 1000

is

n = 110611511026604935641074705584421138393028001852577373936470952377218354575172401275457597579044729873152469512963401398362087144972181770571895264066114088968182356842977823764462179821981744448731785408629116321919957856034605877855212667092287520105386027668843119590555646814038787297694678647529533718769401069269427475868793531944696435696745559289326610132208504257721469829210704462876574915362273129090049477919400226313586033

For this calculation you’ll need to increase the precision from 50 digits to something like 500 digits, something more than 435 because n is a 435-digit number.

In case you’re wondering whether my function for computing harmonic numbers is accurate enough, it’s actually overkill, with error O(1/120n4).

Inverting series that are flat at zero

Posted on by John

The previous post looked at solving the equation

2\theta + \sin(2\theta) = \pi \sin( \psi)

which arises from the Mollweide map projection. Newton’s method works well unless φ is near π/2. Using a modified version of Newton’s method makes the convergence faster when φ = π/2, which is kinda useless because we know the solution is $theta; = π/2 there. When φ is near π/2, the modified Newton’s method may diverge. I ended the previous post by saying a series solution would work better when φ is sufficiently close to π/2. This post will flesh that out.

Let x = π − 2θ. Now the task is to solve

x - \sin(x) = y

for small positive values of y.

The left side is

\frac{x^3}{3!} - \frac{x^5}{5!} + \frac{x^7}{7!} + \cdots

and so for very small values of y, and thus very small values of x, we have

x = (6y)^{1/3}

If this solution is not sufficiently accurate, we can invert the power series above to get a power series in y that gives the solution x. However, the Lagrange inversion theorem does not apply because the series has a zero derivative at 0. Instead, we have to use Puiseux series inversion, looking for a series in y1/3 rather than a series in y. From the Puiseux series we can see that

x = (6y)^{1/3} + y/10

is a more accurate solution. For even more accuracy, you can compute more terms of the Puiseux series.

Fitting a parabola to an ellipse and vice versa

Posted on by John

The previous post discussed fitting an ellipse and a parabola to the same data. Both fit well, but the ellipse fit a little better. This will often be the case because an ellipse has one more degree of freedom than a parabola.

There is one way to fit a parabola to an ellipse at an extreme point: match the end points and the curvatures. This uses up all the degrees of freedom in the parabola.

When you take the analogous approach to fitting an ellipse to a parabola, you have one degree of freedom left over. The curvature depends on a ratio, and so you can adjust the parameters while maintaining the ratio. You can use this freedom to fit the ellipse better over an interval while still matching the curvature at the vertex.

The rest of the post will illustrate the ideas outlined above.

Fitting a parabola to an ellipse

Suppose you have an ellipse with equation

(x/a)² + (y/b)² = 1.

The curvature at (±a, 0) equals a/b² and the curvature at (0, ±b) equals b/a².

Now if you have a parabola

x = cy² + d

then its curvature at y = 0 is 2|c|.

If you want to match the parabola and the ellipse at (a, 0) then d = a.

To match the curvatures at (a, 0) we set a/b² = 2|c|. So c = −a/2b². (Without the negative sign the curvatures would match, but the parabola would turn away from the ellipse.)

Similarly, at (−a, 0) we have d = −a and c = a/2b². And at (0,  ±b) we have d = ±b and c = ∓b/2a².

Here’s an example with a golden ellipse.

Fitting an ellipse to a parabola

Now we fix the parabola, say

y = cx²

and find an ellipse

(x/a)² + ((y − y0)/b)² = 1

to fit at the vertex (0, 0). For the ellipse to touch the parabola at its vertex we must have

((0 − y0)/b)² = 1

and so y0 = b. To match curvature we have

b/2a² = c.

So a and b are not uniquely determined, only the ratio b/a². As long as this ratio stays fixed at 2c, every ellipse will touch at the vertex and match curvature there. But larger values of the parameters will match the parabola more closely over a wider range. In the limit as b → ∞  (keeping b/a²  = 2c), the ellipses become a parabola.

Related posts

Trinomial Coefficients and Kings

Posted on by John

The term trinomial coefficient is used a couple different ways. The most natural use of the term is a generalization of bionomial coefficients to three variables, i.e. numbers of the form

\frac{n!}{i! \, j! \, k!}

where i + j + k = n. These numbers are the coefficients of yj zk in the expansion of (x + y + z)n.

But there’s another use of the term trinomial coefficient, the one that we are interest in for this post, and that is the coefficients of xk in the expansion of (1 + x + 1/x)n. These numbers can be visualized in a triangle analogous to Pascal’s triangle.

In Pascal’s triangle, each entry is the sum of the two entries above it. In the trinomial triangle, each entry is the sum of the three entries above it.

If you’ve been reading this blog regularly you know I’ve been interested in chess puzzles lately. The thing that make me interested in trinomial coefficients is that they relate to moves of a king on a chessboard.

If you note the number paths a king can take to a square using the minimal number of moves, you get the trinomial triangle. The reason is simple: the number of ways to get to a square equals the number of ways to get there via the square up and to the left, plus the number of ways to get their via the square above, plus the number of ways to get their via the square up and to the right.

Related posts

Hand calculating exp(x)

Posted on by John

The previous post mentioned that Martin Gardner announced that Ramanujan’s conjecture that exp(π√163) in an integer had been proven. This was an April Fool’s joke in 1975. Gardner said

Working by hand, he [Ramanujan] found the value to be 262537412640768743.999999999999… The calculations were tedious, and he was unable to verify the next decimal digits.

Calculating exp(π√163) without a computer—Ramanujan died in 1920—would indeed be tedious, but not insurmountable. Certainly it would not stop someone like Ramanujan from testing a conjecture.

How might you go about calculating exp(π√163) by hand?

Algorithm

One possibility is an algorithm in [1].

\exp(t) = \left(1 + r + \frac{r^2}{2!} + \frac{r^3}{3!} + \cdots \right)^{256} 2^n

where r = t′ / 256, t′ = t − n log 2, and n is chosen to minimize |t′|.

We can choose n so that |t′| < log(2)/2 and so |r| < 0.014. This means the infinite series converges rapidly and not too many terms will be needed, depending on the desired precision.

The calculation x256 can be done by squaring 8 times.

Although the context of this post is hand calculations, this would also be a viable algorithm for a program doing extended precision calculations.

Example

In our case, t = π√163 = 40.1091… and we choose n = 58 so that t′ = −0.09336….

Then r = −0.000364….

So each term in the series will contribute 3 or 4 decimal places to the desired precision at first, more once the factorial denominators get large.

Related posts

[1] Jonathan Borwein and David Bailey. Mathematics by Experiment. Volume 1.

Max and min orbital speed

Posted on by John

An earlier post needed to calculate how much the speed of a planet varies in orbit. The planet moves fastest as perihelion, the point in its orbit closest to the sun, and it moves slowest at aphelion, when it is furthest from the sun.

The ratio of the maximum to minimum speed turns out to be a simple expression in terms of the eccentricity e of the orbit:

(1 + e)/(1 − e).

You can derive this fairly quickly from the vis-viva equation, which in turn is derived from conservation of energy.

There are several things I find interesting about this. First, that the expression is so simple. Second, it can be simplified even more for small e:

(1 + e)/(1 − e) ≈ 1 + 2e.

This comes from expanding the ratio as a series:

(1 + e)/(1 − e) = 1 + 2e + 2e² + 2e³ …

This explains two things from the previous post. First, that the variation in orbital speed, for both Earth and Mars, worked out to be about 2e. The eccentricity of Earth’s orbit is 0.0167 and orbital speed varies by about 3%. Mars’ orbit has eccentricity 0.0934 and its orbital speed varies by about 19%. Since the eccentricity of Mars orbit, while fairly small, is larger than that of Earth, the quadratic term matters more for Mars.

Finally, “I keep running into the function f(z) = (1 − z)/(1 + z),” as I first wrote four years ago, and wrote on again a few months ago. It comes up, for example, in computing impedance, in mental calculation tricks, and in efficient calculation of the perimeter of an ellipse. Now you can add to that list calculating the variation in orbital speed in a two body problem.