How to come out with this Overflow in my code?

Question

So, im trying to solve an EDO using Bulirsch-Stoer adaptative method, but im getting an overflower warning. How to solve this?

    from math import sin,cos,pi
from numpy import arange, array, empty, shape,float64
import matplotlib.pyplot as plt
import time

tempo_inicial = time.time()

#Definindo as constantes do exercicio
a = 1
b = 3
ti = 0
tf = 20
prec = 10e-10
x0 = y0 = 0
H = 20
nmax = 8

#Definindo as funções

def f(r,t):
    x = r[0]
    y = r[1]
    f0 = 1-(b+1)*x+a*x**2*y
    f1 = b*x-a*x**2*y
    return(array([f0,f1],float64))

My functions can be seem in f(r,t). I also tried to use numpy float64 but it doesn't change much.

The question in its present form is off-topic, as you ask about the theoretical properties of the ODE system. There is no code presented that could be wrong. The polynomial character of the equation may well lead to positive feed-back and thus to a blow-up of the solution in finite time which would generate overflow errors in any solver. — Lutz Lehmann, May 07 '20 at 09:45
As a dynamical blow-up does not happen when using `scipy.integrate.odeint`, the error has to be in parts of the code that are not shown here. — Lutz Lehmann, May 07 '20 at 09:50
How different is your Bulirsch-Stoer implementation from for instance https://stackoverflow.com/questions/38945920/printing-the-solution-for-bulirsch-stoer-algorithm-in-python? Are you aware that this method is sub-standard in the same way that Euler is inferior to, say, RK4? — Lutz Lehmann, May 09 '20 at 08:23

score 0 · Answer 1 · answered May 09 '20 at 08:13

Not an answer, just expanding on my comment requesting essential details to reconstruct the error: Using scipy.integrate.odeint with the given data and function

from numpy import linspace
from scipy.integrate import odeint

time = linspace(ti, tf, 21)
rsol = odeint(f, [x0,y0], time, atol=prec, rtol=0.01*prec)
for r,t in zip(rsol,time):
    x,y = r;
    print(f't={t:.2f}, x={x:.12f}, y={y:.12f}')

I get without any warning the results

t= 0.00, x=0.000000000000, y=0.000000000000
t= 1.00, x=0.251223552118, y=0.558443392482
t= 2.00, x=0.269384788526, y=1.278992302288
t= 3.00, x=0.285993102849, y=1.984958673784
t= 4.00, x=0.306864172872, y=2.668107280985
t= 5.00, x=0.334793283816, y=3.320108076380
t= 6.00, x=0.375712107368, y=3.925439075135
t= 7.00, x=0.446840813738, y=4.447011243525
t= 8.00, x=0.636678633858, y=4.735903265638
t= 9.00, x=3.584850495692, y=1.405051562362
t=10.00, x=1.764258113375, y=1.397239762022
t=11.00, x=0.648060470900, y=2.400877519264
t=12.00, x=0.384539021730, y=3.190476648398
t=13.00, x=0.379538742491, y=3.821208432875
t=14.00, x=0.436146645718, y=4.361001236766
t=15.00, x=0.588205110586, y=4.711540741504
t=16.00, x=2.339363929033, y=2.958090116701
t=17.00, x=2.049441930465, y=1.253998459089
t=18.00, x=0.751410490789, y=2.250138349553
t=19.00, x=0.397736992172, y=3.081666874339
t=20.00, x=0.375103612719, y=3.727915883209

This looks nicely bounded, no source for your reported error is recognizable.

Note that 10e-10 is 1e-9.

How to come out with this Overflow in my code?

1 Answers1