Forward Euler method for N-dimensional ODE

[Dunnomuch]

Hi,

I'm trying to come up with a FWE method for an N-dimansional ODE. So far, I was able to turn the one dimensional solution into a two dimensional one:

from scitools.std import *
import numpy as numpy
v0=0
K=3 
n=10
dt=0.1
t0=0
 
def f1(u1, u2, u3, t):
    return u1
 
def f2(u1, u2, u3, t):
    return u2
 
def f3(u1, u2, u3, t):
    return u3
 
listf=[f1, f2, f3]
listu0=[1, 1, 1]
listu0 = map(float,listu0)
 
class ForwardEuler:
    def __init__(self, listf, dt):
        self.listf, self.dt = listf, dt
 
    def set_initial_condition(self, listu0, t0=0):
        self.listu = [] 
        self.t = [] 
        self.listu.append(listu0)
        self.t.append(float(t0))
        self.i =0
    def solve(self, K):
        """Advance solution in time until t <= T."""
        tnew = 0
        while tnew <= K:
            listunew = self.advance() 
            self.listu.append(listunew)
            tnew = self.t[-1] + self.dt
            self.t.append(tnew)
            self.i += 1
        return numpy.array(self.listu), numpy.array(self.t)
 
    def advance(self):
        """Advance the solution one time step."""
        # avoid "self." to get more readable formula:
        listu, dt, listf, i, t = \
            self.listu, self.dt, self.listf, self.i, self.t[-1]
 
        lustunew = listu[i] + dt*listf(listu[i], listu[i], t) # Forward Euler scheme
        return listunew
 
method = ForwardEuler(listf, dt)
method.set_initial_condition(listu0, t0)
listf, t = method.solve(K)
print listf
 

I'm now working on the N-dim solution by using a matrix. It looks like this:

class Matrix(object): 
    def __init__(self, kolommen, rijen): 
        self.kolommen = kolommen 
        self.rijen = rijen 
        # Creeer matrix met nullen
        self.matrix = [] 
        for i in range(rijen): 
            ea_rijen = [] 
            for j in range(kolommen): 
                ea_rijen.append(0) 
            self.matrix.append(ea_rijen) 
 
    def plaatselement(self, kolom, rij, v):   #Zet element v op plaats [kolom, rij]
        self.matrix[kolom-1][rij-1] = v 
 
    def neemelement(self, kolom, rij):  #Neemt element op plaats [kolom, rij]
        return self.matrix[kolom-1][rij-1] 
 
    def __repr__(self): #Dient om matrix te printen
        outStr = "" 
        for i in range(self.rijen): 
            outStr += 'Rij %s = %s\n' % (i+1, self.matrix[i]) 
        return outStr 
 

Sorry for the weird names, but I'm from Belgium so I speak dutch and therefor some of the names I use are dutch. Anyway, the beginning of my N-dim solution looks like this:

from scitools.std import *
from numpy import *
from Matrix import Matrix
 
#moeten gedefinieerd worden
F
x0 = []
#N, F, u0 worden gedefinieerd in oproepprogramma
 
class ForwardEuler:
    def __init__(self, F, dt):  
        self.F, self.dt = F, dt
 
    def set_initial_condition(self, Matrix ,lijstU0, N, t0=0):
        M = T/float(self.dt)  #We bepalen de lengte van de kolommen door interval te delen door de lengte van de stapjes
        M = int(M)  #We gebruiken hier int om er een geheel getal van te maken en het rond naar beneden af.
        self.u = Matrix(M,N)  #We maken onze matrix een MxN matrix
 
        n = 1 #Tellen over het aantal functies
        while n <= N:   #Loop plaats cte u0 in eerste kolom
            self.u.plaatselement(1, n, u0[n-1])
            n += 1
 
        self.t = []    
        self.t.append(float(t0))     
        self.i = 0
 

So now I've got to define the functions solve and advance, but I'm a bit stuck here... If anyone could help, that would be great.

[Glenton]

Hi

So numpy has a matrix which you could surely use? It doesn't look like your Matrix class does anything other than assignment and recall?

Could you explain what the n-dimensional step is going to be? In plain english would be fine? What's the calculation that you're struggling with?

I'm relatively confident that numpy matrix or array operations will be suitable for you.

[Dunnomuch]

An N-dimensional ODE has N equasions:

df1/dt = f1(x1, x2, x3, ..., xN, t)

df2/dt = f2(x1, x2, x3, ..., xN, t)

...

dfn/dt = fn(x1, x2, x3, ..., xN, t)

The solution should be N different functions and those should come in the matrix. Afterwards, it should be possible to plot those N functions. But don't get me wrong: We won't really find functions, but functionvalues in a certain interval [0, T]. And it's those values that we plot afterwards.
I hope I solved your question.

[Glenton]

Hi

So this example might get you started (apologies for the scrapiness - I am in haste):

def f0(u,t):
    #u is an array of n elements
    return 2+u[0]+u[1]+t
 
def f1(u,t):
    #u is an array of n elements
    return 1-u[0]-u[2]-t
 
def f2(u,t):
    #u is an array of n elements
    return 3+u[2]+2*u[1]+t
 
 
F=[f0,f1,f2]
U=[0.0,0.0,0.0]
 
results=[]
results.append(U)
 
dt=0.1
t=0
 
T=1
 
while t<=T:
    Utemp=[]
    for u,f in zip(U,F):
        Utemp.append(u+dt*f(U,t))
    t+=dt
    U=Utemp[:]
    results.append(U)
 
print results
 

Of course, if the functions were always linear, then they could be represented by a matrix of coefficients, which might make it all a bit easier, but less general.

It might also be worth checking out the scipy ability to do this kind of thing (which is probably far superior to the above!). See for example here

[Dunnomuch]

I'm a little confused... You don't seem to use the Forward Euler methode...?

[Glenton]

It's U_{i+1} = U_i + dt.f(U,t), right?

Anyway, I was just trying to show you how to code an n-dimensional step forward with n functions and n variables, and a t variable. I was trying to get the FEM based on your first post, but if I haven't, perhaps you'll be able to do so and post back...

Let us know how it goes, and if you have any further questions.