我有一个关于使用 RK4 求解二阶微分方程的问题，考虑了对一阶导数的额外约束。我正在做这里显示的示例，并进行一些修改
θ′′(t( +b θ′(t( + c sin(θ(t(( = 0
附加约束是：
当 θ i θ(i+1(<0，>(i+1(=0.9θ′i，
其中 i 是t 的步骤，i+1 是 i 之后的一步。在现实世界中，当位移方向反转时，其速度降低到90%。
向量上，如果y(t( = (θ(t(，ω(t((，则方程为ẏ= f(t，y(，其中f(t，y( = (y₂(t(， −by₂(t( − cos(y₁(t(((。
在这个问题中，如果我想在ω或 θ′(t( 上添加约束(它们是同一件事(，我应该如何修改代码？ 这是我的代码不起作用。附加条件使 θ′ 不连续。以下"自制"解决方案无法正确更新 θ′。我是Python的新手，这是我的第一篇StackOverflow帖子。非常感谢任何指导。
def rungekutta4(f, y0, t, args=()):
n = len(t)
y = np.zeros((n, len(y0)))
y[0] = y0
for i in range(n - 1):
h = t[i+1] - t[i]
if y[i][0]*y[i+1][0]<0:
k1 = f(y[i], t[i], *args)
k2 = f(y[i] + k1 * h / 2., t[i] + h / 2., *args)
k3 = f(y[i] + k2 * h / 2., t[i] + h / 2., *args)
k4 = f(y[i] + k3 * h, t[i] + h, *args)
y[i+1] = y[i] + (h / 6.) * (k1 + 2*k2 + 2*k3 + k4)*0.9
else:
k1 = f(y[i], t[i], *args)
k2 = f(y[i] + k1 * h / 2., t[i] + h / 2., *args)
k3 = f(y[i] + k2 * h / 2., t[i] + h / 2., *args)
k4 = f(y[i] + k3 * h, t[i] + h, *args)
y[i+1] = y[i] + (h / 6.) * (k1 + 2*k2 + 2*k3 + k4)
return y
除非我完全误解了你，否则你想要的是 f 中的简单大小写区分：在数学上，你有 f₂(t，y( = −by₂(t( −cos(y₁(t(( 如果θi²−1 =y₁(t(²−1 <0 和 0.9·(y₂−1( 否则。这一切仍然只是f对y的依赖，它可以简单地在编程方面实现。
例如，您可以定义：
def f(y):
θ = y[0]
ω = y[1]
return [
θ,
-b*ω-cos(θ) if θ**2<1 else 0.9*(ω-1)
]
虽然这可能按原样工作，但由于f不连续(或具有连续导数(，您可能会遇到问题，特别是如果您想使用具有步长控制的更高级的积分器而不是您自制的积分器。 为避免这种情况，请将if–else构造替换为(尖锐的(sigmoid。有关这方面的更多详细信息，请参阅我的这个答案。
在当前的公式中，假设钟摆每次通过垂直位置时，其速度都会降低 10%，这可以大致安排为
for i in range(n - 1):
h = t[i+1] - t[i]
y[i+1] = RK4step(f,t[i],y[i],h, args)
if y[i+1,0]*y[i,0] < 0: y[i+1,1] *= 0.9
return y
也就是说，首先计算新值，然后应用条件。时间步长应该足够小，角度只能改变几度。对于较大的时间步长，您必须拆分包含零交叉的步骤，使用一些根查找方法(如割线方法(来查找更准确的根时间。