Jul-08-2018, 09:31 PM
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:
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
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))/epsNow 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