在 python 中实现欧拉法求解 ODE

Question

我正在做一个小科学项目，用 ODE（常微分方程）模拟机械系统。

我想使用 Euler 方法来这样做，但是我相信我正在做一些工作，因为我创建的输出数据图不是它应该是的，我是 100%确保等式本身是正确的。你能看看我的代码并告诉我我做错了什么吗？此代码会将 CSV 值输出到控制台和 output.csv 文件中。

测试：

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

对机械系统使用 Euler 通常不是一个好主意。探索此语句的最简单测试用例是简单的振荡器 x''+x=0，您会发现系统的能量快速增长。

对于一般的机械系统，您有一个运动方程 m*x'' = F(t,x,x')。这给你一个矢量值系统

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

注意系统函数中没有h或dt。这对于 ODE 及其数值处理至关重要。轨迹遵循的方向场不依赖于数值求解方法的细节。

因此，请将您的方程重新表述为二阶方程或一阶相位系统 space。

我认为你的等式应该是

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

质量

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

或其某些缩放版本。所以某种 spring 在重力和摩擦力下。如果这是真的，那么您将物理系统的加速度视为数值实现中的速度。这只能给出无意义的结果。

剥离掉类，系统的二阶可以这样反映：

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)

在 python 中实现欧拉法求解 ODE

Implementing Euler's Method in python to solve ODE

python

python-2.7

numerical-methods

numerical

differential-equations