Implementing Euler's Method in python to solve ODE

Question

I am working on a little science project to simulate a mechanical system with ODE (Ordinary Differential Equation).

I would like to use Euler's Method to do so, however I believe that I am doing something work, because the plot that I create of the output data isn't what it is supposed to be and I am 100% sure that the equation itself is correct. Could you look at my code and tell me what I am doing wrong? This code will output CSV values into console and into output.csv file.

Test:

class Calculus:
    def __init__(self):
        print "object created"

    @staticmethod
    def euler(f,y0,a,b,h):
        """y0 - temperatura poczatkowa, a-chwila poczatkowa czasu(warunek brzegowy) , b-chwila koncowa, h - krok"""
        t,y = a,y0
        time = [t]
        value = [y]

        while t <= b:
            print "%6.3f,%6.3f" % (t,y)

            t += h
            y += h * f(t,y, h, value)
            time.append(t)
            value.append(y)

        data = {'time' : time, 'value' : value}
        return data

class Mechanical_system:
    #Constructor
    def __init__(self, momentum1, momentum2, b1, b2, N1, N2, w1, Tau, k):
        self.system_parameters = {'momentum1': momentum1, 'momentum2': momentum2,
                                               'b1': b1, 'b2' : b2, 
                                               'N1': N1, 'N2' : N2, 
                                               'w1' : w1,
                                               'Tau' : Tau,
                                               'k' : k};

    def ODE(self, time, w, h, value):
        """ODE - ordinary differential equation describing our system"""

        #first calculation will only have one value in the list, so we can't calculate delt_value = current - last
        #thus, we need to assume that if there is no penault value (index error) let penault = 0
        try:
            penault = value[len(value) -2]
        except IndexError:
            penault = 0


        momentum1 = self.system_parameters['momentum1']
        momentum2 = self.system_parameters['momentum2']
        b1 =        self.system_parameters['b1']
        b2 =        self.system_parameters['b2']
        N1 =        self.system_parameters['N1']
        N2 =        self.system_parameters['N2']
        Tau =       self.system_parameters['Tau']
        k =         self.system_parameters['k']

        dOmega = w - penault
        dt = h
        Omega = float(dOmega) / float(dt)

        return (1.0 / (momentum1 + ((N1/N2)**2)*momentum2))*(Tau - (Omega)*(b1 + ((N1/N2)**2)*b2) - w * ((N1/N2)**2)*k)

if __name__ == "__main__":
        """reads values from input boxes and then calculates the value and plot value(time)"""
        momentum1 = 1
        momentum2 = 2
        b1 =        0.5
        b2 =        0.6
        N1 =        10
        N2 =        20
        a =         0
        b =         100
        h =         0.1
        w1 =        0
        Tau =       100
        k =         0.5

        system1 = Mechanical_system(momentum1, momentum2, b1, b2, N1, N2, w1, Tau, k)
        data = Calculus.euler(system1.ODE, w1, a, b, h)

        ##writing output to CSV file
        f = open('output.csv', 'w')
        for index in range(len(data['time'])):
            f.write("%s,%s\n" % (data['time'][index], data['value'][index]))

        f.close()

        del system1
        del data

Answer 1

Using Euler for mechanical systems is in general a bad idea. The easiest test case to explore this statement is the simple oscillator x''+x=0 where you will find that the energy of the system grows rapidly.

For a general mechanical system you have an equation of motion m*x'' = F(t,x,x'). This gives you a vector valued system

def f(t,y):
    x,v = y
    return [ v; 1/m*F(t,x,v) ]

Notice that there is no h or dt in the system function. This is essential for ODE and their numerical treatment. The direction field that the trajectories follow do not depend on details of the numerical solution method.

So please reformulate your equation as a second order equation or a first order system in phase space.

---

I take it that your equation is meant to be

F(x,v) = Tau - (b1 + ((N1/N2)**2)*b2) * v -  ((N1/N2)**2)*k * x

with mass

m = momentum1 + ((N1/N2)**2)*momentum2

or some scaled version thereof. So some kind of spring under gravity and friction. If this is true, then you are treating the acceleration of the physical system as the velocity in the numerical implementation. Which can only give meaningless results.

---

Stripping the classes out, the second order of the system could be reflected like this:

def euler(f,y0,v0,a,b,h):
    t,y,v = a,y0,v0
    time = [t]
    value = [y]

    while t <= b:
        print "%6.3f,%6.3f" % (t,y)
        accel =  f(t,y,v)
        t += h
        y += h * v
        v += h*accel
        time.append(t)
        value.append(y)

    data = {'time' : time, 'value' : value}
    return data

def ODE(t,y,v,param):
    momentum1, momentum2, b1, b2, N1, N2, Tau, k = param
    mass = momentum1 + ((N1/N2)**2)*momentum2
    force = Tau - v*(b1 + ((N1/N2)**2)*b2) - y * ((N1/N2)**2)*k
    return force/mass


momentum1 = 1
momentum2 = 2
b1 =        0.5
b2 =        0.6
N1 =        10.
N2 =        20.
Tau =       100
k =         0.5    
a =         0
b =         100
h =         0.1
y0 =        0
v0 =        0

params = [momentum1, momentum2, b1, b2, N1, N2, Tau, k]
data = euler(lambda t,y,v: ODE(t,y,v,params),y0,v0,a,b,h)