用四阶龙格库塔求解阻尼振荡器的二阶微分方程
Using 4th order Runge Kutta to solve the 2nd order differential equation of a damped oscillator
所以我的任务是使用四阶龙格库塔方法求解阻尼振荡器的二阶微分方程。
我的 runge-kutta 方法的函数看起来是这样的
def RungeKutta(f,y0,x):
y=np.zeros((len(x),len(y0)))
y[0,:]=np.array(y0)
h=x[1]-x[0]
for i in range(0,len(x)-1):
k1=h*np.array(f(y[i,:],x[i]))
k2=h*np.array(f(y[i,:]+k1/2,x[i]+h/2))
k3=h*np.array(f(y[i,:]+k2/2,x[i]+h/2))
k4=h*np.array(f(y[i,:]+k3,x[i]+h))
y[i+1,:]=y[i,:]+k1/6+k2/3+k3/3+k4/6
return y
rungeKutta 函数工作正常，我已经用示例输入列表对其进行了测试，所以这似乎不是问题所在
我给了
question parameters
并且必须做一个class来解决问题
class harmonicOscillator:
def __init__(self,m,c,k):
if((m>0) and ((type(m) == int) or (type(m) == float))):
self.m = m
else:
raise ValueError
if((c>0) and ((type(c) == int) or (type(c) == float))):
self.c = c
else:
raise ValueError
if((k>0) and ((type(k) == int) or (type(k) == float))):
self.k = k
else:
raise ValueError
def period(self):
self.T = 2 * np.pi * (self.m / self.k)**(0.5)
return(self.T)
def solve(self, func, y0):
m = self.m
c = self.c
k = self.k
T = self.T
t = np.linspace(0,10*T,1000)
但我不确定真正的进步在哪里。我试过将二阶微分方程变成这样的 lambda 函数
F = lambda X,t: [X[1], (-c) * X[1] + (-k) * X[0] + func(t)]
然后将其传递到我的 RungeKutta 函数中
result = RungeKutta(F, y0, t, func)
return(result)
但我不是很精通 lambda 函数，而且显然在某处出错了。
它应该传递的示例输入是这样的
####### example inputs #######
m=1
c=0.5
k=2
a = harmonicOscillator(m,c,k)
a.period()
x0 = [0,0]
tHO,xHO= a.solve(lambda t: omega0**2,x0)
非常感谢您的帮助。问题的要求是我必须使用上面的 rungeKutta 函数，但我现在有点迷路了
谢谢。
我认为外部强迫项和 Runge Kutta 导数辅助函数可能存在一些混淆 F。 RK4中的Freturns一阶微分方程组X的导数dX/dt。不幸的是，阻尼振荡器中的强迫项也称为 F，但它是 t.
的函数
您的一个问题是 RungeKutta() 函数的元数（参数数量）与对该函数的调用不匹配：您尝试执行 RungeKutta(F, y0, t, func)，但 RungeKutta() 函数只接受按顺序 (f, y0, x) 的参数。
换句话说，你当前RungeKutta()函数中的f参数应该封装强制函数F(t)。
你可以通过助手来做到这一点：
# A constant function in your case, but this can be any function of `t`
def applied_force(t):
# Note, you did not provide a value for `omega0`
return omega0 ** 2
def rk_derivative_factory(c, k, F):
return lambda X, t: np.array([X[1], -c * X[1] - k * X[0] + F(t)])
rk_derivative_factory() 是一个函数，它采用阻尼系数、spring 常数和一个强制函数 F(t)，returns function 以系统 X 和时间步长 t 作为参数（因为这是 RungeKutta() 的实现所要求的）。
然后，
omega0 = 0.234
m, c, k = 1, 0.25, 2
oscillator = HarmonicOscillator(m, c, k)
f = rk_derivative_factory(oscillator, applied_force)
x_osc = oscillator.solve(f, [1, 0])
其中 solve() 定义如下：
def solve(self, func, y0):
t = np.linspace(0,10 * self.period(), 1000)
return RungeKutta(f, y0, t)
顺便说一句，我强烈建议您的变量名称更加一致。您将振荡器的初始状态命名为 x0，并将其作为参数 y0 的参数传递给 RungeKutta()，以及 RungeKutta() 的 x 参数代表时间...变得很混乱。
完整解决方案
最后，您对 RK4 的实现不太正确，所以我已经解决了这个问题并做了一些其他的小改进。
请注意，您可能要考虑的一件事是让 HarmonicOscillator.solve() 函数采用求解器。然后你可以尝试不同的集成商。
import numpy as np
def RungeKutta(f, y0, x):
y = np.zeros((len(x), len(y0)))
y[0, :] = np.array(y0)
h = x[1] - x[0]
for i in range(0, len(x) - 1):
# Many slight changes below
k1 = np.array(f(y[i, :], x[i]))
k2 = np.array(f(y[i, :] + h * k1 / 2, x[i] + h / 2))
k3 = np.array(f(y[i, :] + h * k2 / 2, x[i] + h / 2))
k4 = np.array(f(y[i, :] + h * k3, x[i] + h))
y[i + 1, :] = y[i, :] + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
return y
# A constant function in your case, but this can be any function of `t`
def applied_force(t):
# Note, you did not provide a value for `omega0`
return omega0 ** 2
def rk_derivative_factory(osc, F):
return lambda X, t: np.array([X[1], (F(t) - osc.c * X[1] - osc.k * X[0]) / osc.m])
class HarmonicOscillator:
def __init__(self, m, c, k):
if (type(m) in (int, float)) and (m > 0):
self.m = m
else:
raise ValueError("Parameter 'm' must be a positive number")
if (type(c) in (int, float)) and (c > 0):
self.c = c
else:
raise ValueError("Parameter 'c' must be a positive number")
if (type(k) in (int, float)) and (k > 0):
self.k = k
else:
raise ValueError("Parameter 'k' must be a positive number")
self.T = 2 * np.pi * (self.m / self.k)**(0.5)
def period(self):
return self.T
def solve(self, func, y0):
t = np.linspace(0, 10 * self.period(), 1000)
return RungeKutta(func, y0, t)
演示：
import plotly.graph_objects as go
omega0 = 0.234
m, c, k = 1, 0.25, 2
oscillator = HarmonicOscillator(m, c, k)
f = rk_derivative_factory(oscillator, applied_force)
x_osc = oscillator.solve(f, [1, 0])
x, dx = x_osc.T
t = np.linspace(0, 10 * oscillator.period(), 1000)
fig = go.Figure(go.Scatter(x=t, y=x, name="x(t)"))
fig.add_trace(go.Scatter(x=t, y=dx, name="x'(t)"))
输出：
所以我的任务是使用四阶龙格库塔方法求解阻尼振荡器的二阶微分方程。
我的 runge-kutta 方法的函数看起来是这样的
def RungeKutta(f,y0,x):
y=np.zeros((len(x),len(y0)))
y[0,:]=np.array(y0)
h=x[1]-x[0]
for i in range(0,len(x)-1):
k1=h*np.array(f(y[i,:],x[i]))
k2=h*np.array(f(y[i,:]+k1/2,x[i]+h/2))
k3=h*np.array(f(y[i,:]+k2/2,x[i]+h/2))
k4=h*np.array(f(y[i,:]+k3,x[i]+h))
y[i+1,:]=y[i,:]+k1/6+k2/3+k3/3+k4/6
return y
rungeKutta 函数工作正常，我已经用示例输入列表对其进行了测试，所以这似乎不是问题所在
我给了 question parameters 并且必须做一个class来解决问题
class harmonicOscillator:
def __init__(self,m,c,k):
if((m>0) and ((type(m) == int) or (type(m) == float))):
self.m = m
else:
raise ValueError
if((c>0) and ((type(c) == int) or (type(c) == float))):
self.c = c
else:
raise ValueError
if((k>0) and ((type(k) == int) or (type(k) == float))):
self.k = k
else:
raise ValueError
def period(self):
self.T = 2 * np.pi * (self.m / self.k)**(0.5)
return(self.T)
def solve(self, func, y0):
m = self.m
c = self.c
k = self.k
T = self.T
t = np.linspace(0,10*T,1000)
但我不确定真正的进步在哪里。我试过将二阶微分方程变成这样的 lambda 函数
F = lambda X,t: [X[1], (-c) * X[1] + (-k) * X[0] + func(t)]
然后将其传递到我的 RungeKutta 函数中
result = RungeKutta(F, y0, t, func)
return(result)
但我不是很精通 lambda 函数，而且显然在某处出错了。
它应该传递的示例输入是这样的
####### example inputs #######
m=1
c=0.5
k=2
a = harmonicOscillator(m,c,k)
a.period()
x0 = [0,0]
tHO,xHO= a.solve(lambda t: omega0**2,x0)
非常感谢您的帮助。问题的要求是我必须使用上面的 rungeKutta 函数，但我现在有点迷路了 谢谢。
我认为外部强迫项和 Runge Kutta 导数辅助函数可能存在一些混淆 F。 RK4中的Freturns一阶微分方程组X的导数dX/dt。不幸的是，阻尼振荡器中的强迫项也称为 F，但它是 t.
您的一个问题是 RungeKutta() 函数的元数（参数数量）与对该函数的调用不匹配：您尝试执行 RungeKutta(F, y0, t, func)，但 RungeKutta() 函数只接受按顺序 (f, y0, x) 的参数。
换句话说，你当前RungeKutta()函数中的f参数应该封装强制函数F(t)。
你可以通过助手来做到这一点：
# A constant function in your case, but this can be any function of `t`
def applied_force(t):
# Note, you did not provide a value for `omega0`
return omega0 ** 2
def rk_derivative_factory(c, k, F):
return lambda X, t: np.array([X[1], -c * X[1] - k * X[0] + F(t)])
rk_derivative_factory() 是一个函数，它采用阻尼系数、spring 常数和一个强制函数 F(t)，returns function 以系统 X 和时间步长 t 作为参数（因为这是 RungeKutta() 的实现所要求的）。
然后，
omega0 = 0.234
m, c, k = 1, 0.25, 2
oscillator = HarmonicOscillator(m, c, k)
f = rk_derivative_factory(oscillator, applied_force)
x_osc = oscillator.solve(f, [1, 0])
其中 solve() 定义如下：
def solve(self, func, y0):
t = np.linspace(0,10 * self.period(), 1000)
return RungeKutta(f, y0, t)
顺便说一句，我强烈建议您的变量名称更加一致。您将振荡器的初始状态命名为 x0，并将其作为参数 y0 的参数传递给 RungeKutta()，以及 RungeKutta() 的 x 参数代表时间...变得很混乱。
完整解决方案
最后，您对 RK4 的实现不太正确，所以我已经解决了这个问题并做了一些其他的小改进。
请注意，您可能要考虑的一件事是让 HarmonicOscillator.solve() 函数采用求解器。然后你可以尝试不同的集成商。
import numpy as np
def RungeKutta(f, y0, x):
y = np.zeros((len(x), len(y0)))
y[0, :] = np.array(y0)
h = x[1] - x[0]
for i in range(0, len(x) - 1):
# Many slight changes below
k1 = np.array(f(y[i, :], x[i]))
k2 = np.array(f(y[i, :] + h * k1 / 2, x[i] + h / 2))
k3 = np.array(f(y[i, :] + h * k2 / 2, x[i] + h / 2))
k4 = np.array(f(y[i, :] + h * k3, x[i] + h))
y[i + 1, :] = y[i, :] + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
return y
# A constant function in your case, but this can be any function of `t`
def applied_force(t):
# Note, you did not provide a value for `omega0`
return omega0 ** 2
def rk_derivative_factory(osc, F):
return lambda X, t: np.array([X[1], (F(t) - osc.c * X[1] - osc.k * X[0]) / osc.m])
class HarmonicOscillator:
def __init__(self, m, c, k):
if (type(m) in (int, float)) and (m > 0):
self.m = m
else:
raise ValueError("Parameter 'm' must be a positive number")
if (type(c) in (int, float)) and (c > 0):
self.c = c
else:
raise ValueError("Parameter 'c' must be a positive number")
if (type(k) in (int, float)) and (k > 0):
self.k = k
else:
raise ValueError("Parameter 'k' must be a positive number")
self.T = 2 * np.pi * (self.m / self.k)**(0.5)
def period(self):
return self.T
def solve(self, func, y0):
t = np.linspace(0, 10 * self.period(), 1000)
return RungeKutta(func, y0, t)
演示：
import plotly.graph_objects as go
omega0 = 0.234
m, c, k = 1, 0.25, 2
oscillator = HarmonicOscillator(m, c, k)
f = rk_derivative_factory(oscillator, applied_force)
x_osc = oscillator.solve(f, [1, 0])
x, dx = x_osc.T
t = np.linspace(0, 10 * oscillator.period(), 1000)
fig = go.Figure(go.Scatter(x=t, y=x, name="x(t)"))
fig.add_trace(go.Scatter(x=t, y=dx, name="x'(t)"))
输出：