Euler class for solve System of ODE

[drudox]

Hi all

I've defined a python class in order to compute the solution of system of differential eq. Do ding so I define a classes named Rhs (right and side) that should represent the right and side of the dy/dt(i-th) this class contains a single float value (initial time , initial value , final time ) and a function (array of function) in order to define this array I simply defined 3 lambda function that rapresent equation(i) and create a np.array of this function

here the Rhs class:

 import types
import numpy as np
import matplotlib.pyplot as plt 


class Rhs:
   '''
      class to define a differential problems 
      contains the dy/dt function and the initial value 
      in order to close the ode problem
   '''
   solution = None

   def __init__(self, fnum : np.ndarray , t0: np.float, tf: np.float, y0 : np.array, n: int , fanal):
      '''   
         Input : 
            - fnum : Function f(t,y(t)) = dy/dt (numeric)
            - 
            - Initial Time t0 
            - Final Time tf
            - Initial value y0(t0) 
      '''
      self.func = fnum
      Rhs.solution  = fanal
      self.t0   = t0
      self.tf   = tf
      self.n    = n
      self.u0   = y0


   def createArray(self):
      '''   
         Create the Array time and f(time) 
         - the time array can be create now 
         - the solution array can be just initialize with the IV
      ''' 
      self.t = np.linspace(self.t0, self.tf, self.n )
      self.u = np.array([self.u0  for i in range(self.n) ])
      
      
      return self.t,self.u   
   
   def f(self,ti,ui):
       return  np.array([function(ti,ui) for function in self.func])     
       

   def Df(self,ti,ui):
      eps = 10e-6
      return ((self.func(ti,ui)+eps) - self.f(ti,ui))/eps

Now I've defined the Euler class :

class Euler(metaclass=ABCMeta):
   ''' 
      Class base for the two Euler solvers (Explicit and Implicit)
   '''
   def __init__(self, dydt : rhs.Rhs) :
      self.dydt   = dydt
      self.dt     = (dydt.tf-dydt.t0)/dydt.n
   
   def plot(self,*args):
       '''  
         Plot the result of the sub classes integrations
         Implicit and Explicit method         
       ''' 
       self.parameter= [*args]
       plt.style.use(['mystyle', 'mystyle-nb'])
       ax = plt.subplot(111)
       ax.plot(self.time,self.u)
       ax.set_title (self.parameter[0])
       ax.set_xlabel(self.parameter[1])
       ax.set_ylabel(self.parameter[2])
       plt.show()

   @abstractmethod
   def solve(self):    
      pass 
   
#---------------------------------------------------------------------------------------------
#     
#      
#
#
#---------------------------------------------------------------------------------------------

class Explicit(Euler):
   '''
      Solve the ODE using first order Foward difference
      Explicit first order Method (Euler)  
   '''
   
   solved = False  # True if the method 'solve' was called without error 

   def __init__(self, dydt: rhs.Rhs, save : bool=True, _print: bool=False,  filename : str=None):
      ''' 
         Initialize the Foward Euler solver 
            - dydt (to the super class) : RHS problem dy/dt = f(t,y(t)) t(0)=y0 
            - save : if True returns the 2 vector time and u = du/dt 
      '''
      self._print = _print
      self.save   = save 
      self.file   = filename
    
      super().__init__(dydt) 

   def solve(self):
      self.time, self.u = self.dydt.createArray() 
      for i in range(len(self.time)-1):
          #for j in range(len(self.u[0,:])):
              self.u[i+1] = self.u[i] + self.dt*self.dydt.f(self.time[i],self.u[i])
         
          
      Explicit.solved = True 
      
      if self._print:
          with open(self.file,'w') as f:
              for i in range(len(self.time)):      
                 f.write('%.4f %4f %4f %4f\n' %(self.time[i] ,self.u[i,0], self.u[i,1], self.u[i,2]))
      
      if self.save:
         return self.time,self.u

   def plot(self): 
     '''
         Check if the problem is solved (if the solve method was already called) 
         and if it is plot the numerical solution of the IV differential problem
     '''
     if Explicit.solved:
         super().plot('ODE Solution using Foward Euler', 'time [s]', 'y(t)')  
     else:    
         print("Unsolved problem, call `solve` method before")

Now in order test the class I've wrote the Lorentz attractor problem in a main file:

import numpy as np
import rhs
import euler

## Lorentz attractor system

func1 = lambda t,u : 10   * (u[1] - u[0])
func2 = lambda t,u : 28   * u[0] - u[1] - u[0] * u[2]  
func3 = lambda t,u : -8/3 * u[2] + u[0]*u[1]

#-------------------------------------------------------------
func = np.array([func1,func2,func3])
#

def main():

   y0 = np.array([1.,0.,0.]) 
 
   problem = rhs.Rhs(func,0.0,100.0,y0,1000) 
   t,u       = problem.createArray()
   fwdeuler_p1 = euler.Explicit(problem3, True,_print=True, filename ='lorentz.dat')
   
   fet,feu = fwdeuler_p1.solve()

if __name__ == '__main__':
    main()

Now, if I try to run this I got this error:

drive.py:31: RuntimeWarning: overflow encountered in double_scalars
  func2 = lambda t,u : 28   * u[0] - u[1] - u[0] * u[2]
drive.py:32: RuntimeWarning: overflow encountered in double_scalars
  func3 = lambda t,u : -8/3 * u[2] + u[0]*u[1]
drive.py:31: RuntimeWarning: invalid value encountered in double_scalars
  func2 = lambda t,u : 28   * u[0] - u[1] - u[0] * u[2]
/home/marco/Programming/Python/Numeric/OdeSystem/euler.py:77: RuntimeWarning: invalid value encountered in add
  self.u[i+1] = self.u[i] + self.dt*self.dydt.f(self.time[i],self.u[i])

It is very strange that I got this problem because this is the same logic that I used developing a fine working code in C++ (if you are interesting to have a look you can go here and find the c++ hierarchy )

I hope that the logic it's clear basically using (in this case) a 3 * 1000 matrix u1,u2,u3 solution vectors (once for each equations) I'm passing all the three vectors to the farray function in this way is like if the 3 equations are coupling and for each u[i] correspond one vector ui ... I hope that somebody help me ... I'm quite noobs in python I have experience with C++ , C and Fortran.
thanks in advance