Making the Madhava series practically useful

The Madhava series is slow and boring, says … Madhava

यत्सङ्ख्ययात्र हरणे कृते निवृत्ता हृतिस्तु जामितया
तस्या ऊर्ध्वगता या समसङ्ख्या तद्दलं गुणोऽन्ते स्यात् ||
तद्वर्गो रूपयुतो हारो व्यासाब्धिघातात् प्राग्वत् 
ताभ्यामाप्तं स्वमृणे कृते धने क्षेप एव करणीयः ||
लभ्धः परिधिः सूक्ष्मो बहुकृत्वो हरणतोऽतिसूक्ष्मश्च ||
When you stop division out of boredom, remember the last divisor
take the next even number, and halve it. That is your numerator.
The square of the numerator, plus one is your divisor. Multiply by four times the diameter as before.
The correction you thus calculate, add to the result if the last operation you did was subtraction, and vice versa.
The corrected diameter you get is accurate, much more so that by a lot of division (ie: by taking several more terms!)

What exactly does this mean?

We are back with व्यासे वारिधिनिहते, the Madhava Pi (or more accurately circumference) series. We saw in an earlier article how it was demonstrated to be correct by Yuktibhāṣā. To jog the reader's memory, the verse is :

व्यासे वारिधिनिहते रूपहृते व्याससागराभिहते ।
त्रिशरादिविषमसंख्याभक्तमृणं स्वं पृथक्क्रमात् कुर्यात् ॥

which is equivalent to

\( C = 4.d – \frac{4.d}{3} + \frac{4.d}{5} – \frac{4.d}{7} + … \)
\( \frac{\pi}{4} = 1 – \frac{1}{3} + \frac{1}{5} – \frac{1}{7} \)

The one problem with this series is that it converges extremely slowly. This is precisely why the verse says निवृत्ता हृतिस्तु जामितया (“when you stop dividing out of boredom”). Remember that all this calculation is happening by hand. How many terms would we calculate? 10? 20? 50?

3.14159265358979

These are the first 15 digits of Pi – this is our target.

Madhava's series converges really really slowly towards this number. Setting the diameter to 1, we add more and more terms to the series to calculate the circumference as

TermsResult
103.0418396189294024
203.0916238066678385
303.1082685666989462
403.1165965567938323
503.1215946525910105
1003.131592903558553
Slow convergence of the Madhava series

50 terms is a lot to calculate by hand, and we haven't got anywhere close to the answer we want. Even after a 100 terms, Aryabhata's 3.1416 is still much better. This series as is, is a theoretical curiosity, and is not much use for practical calculation. Or is it?
(We have conveniently ignored that Madhava himself would have calculated in fractions, not decimals, and would've probably set d = 100000000000, but that doesn't change our point)

TermsResultError1/Terms
103.04183961892940240.099753034660390760.1
203.09162380666783850.049968846921954580.05
303.10826856669894620.033324086890846870.03333..
403.11659655679383230.024996096795960820.025
503.12159465259101050.0199980009987825720.02
1003.1315929035585530.0099997500312403180.01
Error Pattern

We notice that the difference between Pi and our approximation is almost the same, but a little less than (and grows increasingly close to) the reciprocal of the number of terms we have calculated. If only we could add a correction term …

संस्कारः – The correction term.

Let's start with the presumption that we want a correction term, to be added to the end of the series, depending on the number of terms in the series. The correction term must

  1. Be dependent on the number of terms already calculated (not a constant!)
  2. Approximate the rest of the series
  3. Converge faster than the series itself.

Let's assume we have computed n terms in the series. The last denominator we computed was \(p = 2n-1\). We are looking for a correction term of the form \(1/a_p\) to add to the series, such that
\( C \approx 4.d – \frac{4.d}{3} + \frac{4.d}{5} – \frac{4.d}{7} + … (-1)^{\frac{p-1}{2}}.\frac{4.d}{p} + (-1)^{\frac{p+1}{2}}.\frac{4.d}{a_p} \)
As n (and p) grow larger, we want the series to converge faster with the correction term added. ie, we want:
\( 4.d – \frac{4.d}{3} + \frac{4.d}{5} – \frac{4.d}{7} + … (-1)^{\frac{p-3}{2}}.\frac{4.d}{p-2} + (-1)^{\frac{p-1}{2}}.\frac{4.d}{a_{p-1}} \approx 4.d – \frac{4.d}{3} + \frac{4.d}{5} – \frac{4.d}{7} + … (-1)^{\frac{p-1}{2}}.\frac{4.d}{p} + (-1)^{\frac{p+1}{2}}.\frac{4.d}{a_p} \) which, after accounting for signs, simplifies to:
\( \frac{1}{a_{p-1}} + \frac{1}{a_{p}} \approx \frac{1}{p}\)

We define the स्थौल्य (error term) \( E_p = \frac{1}{a_{p-1}} + \frac{1}{a_{p}} – \frac{1}{p} \), which we want to keep as small as we can, and want it to go to zero with n much faster than \(\frac{1}{n}\) does.

First order correction term

Our first attempt at a correction term is \( a_p = 2(p+1) = 4n \). This is easy to calculate. Whatever your last denominator was, add one to it, and double your result to get the denominator of your correction term. This is a good start, as we can see :

TermsResult with first order correction
103.141839618929402
203.141623806667839
303.1416019000322795
403.141596556793832
503.1415946525910106
1003.1415929035585526
Improved convergence of the Madhava series

We are as good as Aryabhata after 20 terms. In 40 terms, we are ahead of him by one digit. By term 100, we have one more digit of accuracy. Good, but we need to do better.

We note that the स्थौल्य is now
\( \begin{align*}
E_p &= \frac{1}{a_{p-1}} + \frac{1}{a_{p}} – \frac{1}{p} \\
&= \frac{1}{2(p-1)} + \frac{1}{2(p+1)} – \frac{1}{p} \\
&= \frac{1}{p^3-p}
\end{align*}\)

We note that a positive स्थौल्य indicates an over-correction, which we can verify from the table. It goes down as the third power of p, indicating that we now converge much faster with this correction term, as we can see from the table.

Second order correction

To improve the correction term, we therefore have to increase its denominator. Can we increase it by a constant 1? That will give us

\( \begin{align*}
E_p &= \frac{1}{a_{p-1}} + \frac{1}{a_{p}} – \frac{1}{p} \
&= \frac{1}{2p-1} + \frac{1}{2p+3} – \frac{1}{p} \
&= \frac{-2p+3}{4p^3+4p^2-3p}
\end{align*}\)

We have now ended up with a term proportional to p in the error numerator. Therefore, Yuktibhāṣā reasons that the denominator of the end-correction is better increased by a factor of 1/p rather than 1, and arrives at the optimal value as \( a_p = 2p+2 + \frac{4}{2p+2}\) That gives us a स्थौल्य of

\( \begin{align*}
E_p &= \frac{1}{a_{p-1}} + \frac{1}{a_{p}} – \frac{1}{p} \
&= \frac{1}{2p-2+\frac{4}{2p-2}} + \frac{1}{2p+2+\frac{4}{2p+2}} – \frac{1}{p} \
&= \frac{-4}{p^5+4p}
\end{align*}\)

This is a good स्थौल्य – negative, implying a slight under-correction, and reducing as the fifth power of p, which means we converge fast with this correction term. Let's see if that really happens

TermsResult with second order correctionAccuracy (decimal places)
103.1415902423707995
203.1415925761868897
303.14159264334432247
403.14159265115408868
503.14159265279099039
1003.141592653564802411
Better improved convergence of the Madhava series

We have nine decimal places of accuracy with 50 terms, and 11 with 100, which is quite impressive. So, our correction term corresponds to
\( a_p = 2p+2 + \frac{4}{2p+2}\)
The term itself will be proportional to

\( \frac{1}{a_p} = \frac{1}{2p+2 + \frac{4}{2p+2}} = \frac{\frac{p+1}{2}}{(p+1)^2+1}\)

(ie: तस्या ऊर्ध्वगता या समसङ्ख्या तद्दलं गुणोऽन्ते स्यात् || तद्वर्गो रूपयुतो हारो …). Adding this correction term, our series now becomes:

\( C \approx 4.d – \frac{4.d}{3} + \frac{4.d}{5} – \frac{4.d}{7} + … (-1)^{\frac{p-1}{2}}.\frac{4.d}{p} + (-1)^{\frac{p+1}{2}}.\frac{4.d.\frac{p+1}{2}}{(p+1)^2+1} \)

The correction term has the opposite sign to the last term in the series – we add if we have subtracted, and vice versa (स्वमृणे कृते धने क्षेप एव करणीयः).

This form, as the table above attests, is good enough for calculation (लभ्धः परिधिः सूक्ष्मो), already to a degree of accuracy that was not known anywhere in Madhava's time. Comparing this table with the first, we see that it is a lot more useful to compute this correction term than to compute many many more terms of the series (बहुकृत्वो हरणतोऽतिसूक्ष्मश्च).

By improving the end correction term further, we can do even better than this, but we will leave that for a later article.

Leave a comment

Leave a Reply

Close Bitnami banner
Bitnami