Python - Error 34, Result Too large

Question

I’m trying to get a sequence of a recursive function, but it always returns me the overflow error.

import math

f = 1.5;
a = 100;
b = 10;
d = 1;
w = 3.76;

x1 = 0;
y1 = 0;
z1 = 0;

h = 10^-3;

for i in range (20000):
    x1 = (1/6)*h**3*y1*b-(1/2)*h**3*y1*a*x1**2-(1/6)*h**3*d*b*x1+(1/6)*h**3*d**2*y1+(1/6)*h**3*d*a*x1**3-(1/6)*h**3*d*f*math.cos(z1)-(1/6)*h**3*f*math.sin(z1)*w+(1/2)*h**2*b*x1-(1/2)*h**2*d*y1;-(1/2)*h**2*a*x1**3+(1/2)*h**2*f*math.cos(z1)+h*y1+x1;
    y1 = y1+h**3*y1*d*a*x1**2-(1/2)*h**3*a*x1**2*f*math.cos(z1)+(1/6)*h**3*w*d*f*math.sin(z1)+(1/6)*h**3*b**2*x1+(1/2)*h**3*a**2*x1**5-(1/6)*h**3*d**3*y1+(1/2)*h**2*y1*b+(1/2)*h**2*d**2*y1+h*b*x1-h*d*y1-h*a*x1**3+h*f*math.cos(z1)-h**3*a*x1*y1**2;-(1/3)*h**3*y1*d*b-(2/3)*h**3*b*a*x1**3+(1/6)*h**3*d**2*b*x1-(1/6)*h**3*d**2*a*x1**3-(3/2)*h**2*y1*a*x1**2;-(1/2)*h**2*d*b*x1+(1/2)*h**2*d*a*x1**3+(1/6)*h**3*b*f*math.cos(z1)+(1/6)*h**3*d**2*f*math.cos(z1)-(1/6)*h**3*f*math.cos(z1)*w**2;-(1/2)*h**2*d*f*math.cos(z1)-(1/2)*h**2*f*math.sin(z1)*w;
    z1 = h*w+z1;

Answer 1

By default, Python uses 64bit floating point numbers, which are the natives of most modern Cpus.

These numbers are very flexible, but they have their limits - and if any operation breaks for a very large number, or very small, you will get a mistake.

The language also offers the module decimal - which allows decimal numbers with arbitrary precision (operations are done in software) - you can configure decimal. Decimal to provide the precision you need (there is a context where you select the number of decimal places, for example).

However, the decimal library does not have the trigonometric functions, which you are using - it is necessary to use a third library for large numbers then.

A good candidate is "mpmath" - (can be installed with pip install mpmath) - for me generated easy cos(0.5) with 500 decimal places:

In [1]: import mpmath                                                                    

In [2]: mpmath.mp.dps = 500                                                              

In [3]: mpmath.cos(0.5)                                                                  
Out[3]: mpf('0.87758256189037271611628158260382965199164519710974405299761086831595076327421394740579418408468225835547840059310905399341382797683328026679975612095022401558762915687859072347693931098961673967701440899764912857021346821838454381839331616880754066081115940348983190805262434229367983882103953443260971069339648047544648581904315236807834735418729899796204210738598702695348232436617950049022784143745936541443062389248031244110265165332107358580108014739732581731402408144610066382413721269523135836227')