在 python 中实现欧拉法求解 ODE
Implementing Euler's Method in python to solve ODE
我正在做一个小科学项目,用 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
对机械系统使用 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)
我正在做一个小科学项目,用 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
对机械系统使用 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)