我在这个 Dopri5 实现中做错了什么
What am I doing wrong in this Dopri5 implementation
我完全不熟悉python,并尝试整合以下颂歌:
$\dot{x} = -2x-y^2$
$\dot{y} = -y-x^2
这会产生一个数组,但所有内容均为 0
我究竟做错了什么?它主要是复制代码,并且与另一个未耦合的代码一起工作正常。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode
def fun(t, z):
"""
Right hand side of the differential equations
dx/dt = -omega * y
dy/dt = omega * x
"""
x, y = z
f = [-2*x-y**2, -y-x**2]
return f
# Create an `ode` instance to solve the system of differential
# equations defined by `fun`, and set the solver method to 'dop853'.
solver = ode(fun)
solver.set_integrator('dopri5')
# Set the initial value z(0) = z0.
t0 = 0.0
z0 = [0, 0]
solver.set_initial_value(z0, t0)
# Create the array `t` of time values at which to compute
# the solution, and create an array to hold the solution.
# Put the initial value in the solution array.
t1 = 2.5
N = 75
t = np.linspace(t0, t1, N)
sol = np.empty((N, 2))
sol[0] = z0
# Repeatedly call the `integrate` method to advance the
# solution to time t[k], and save the solution in sol[k].
k = 1
while solver.successful() and solver.t < t1:
solver.integrate(t[k])
sol[k] = solver.y
k += 1
# Plot the solution...
plt.plot(t, sol[:,0], label='x')
plt.plot(t, sol[:,1], label='y')
plt.xlabel('t')
plt.grid(True)
plt.legend()
plt.show()
您的初始状态 (z0
) 是 [0,0]
。此初始状态的时间导数 (fun
) 也是 [0,0]
。因此,对于这个初始条件,[0,0]
是所有时间的正确解。
如果将初始条件更改为其他值,您应该会观察到更有趣的结果。
我完全不熟悉python,并尝试整合以下颂歌:
$\dot{x} = -2x-y^2$
$\dot{y} = -y-x^2
这会产生一个数组,但所有内容均为 0 我究竟做错了什么?它主要是复制代码,并且与另一个未耦合的代码一起工作正常。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode
def fun(t, z):
"""
Right hand side of the differential equations
dx/dt = -omega * y
dy/dt = omega * x
"""
x, y = z
f = [-2*x-y**2, -y-x**2]
return f
# Create an `ode` instance to solve the system of differential
# equations defined by `fun`, and set the solver method to 'dop853'.
solver = ode(fun)
solver.set_integrator('dopri5')
# Set the initial value z(0) = z0.
t0 = 0.0
z0 = [0, 0]
solver.set_initial_value(z0, t0)
# Create the array `t` of time values at which to compute
# the solution, and create an array to hold the solution.
# Put the initial value in the solution array.
t1 = 2.5
N = 75
t = np.linspace(t0, t1, N)
sol = np.empty((N, 2))
sol[0] = z0
# Repeatedly call the `integrate` method to advance the
# solution to time t[k], and save the solution in sol[k].
k = 1
while solver.successful() and solver.t < t1:
solver.integrate(t[k])
sol[k] = solver.y
k += 1
# Plot the solution...
plt.plot(t, sol[:,0], label='x')
plt.plot(t, sol[:,1], label='y')
plt.xlabel('t')
plt.grid(True)
plt.legend()
plt.show()
您的初始状态 (z0
) 是 [0,0]
。此初始状态的时间导数 (fun
) 也是 [0,0]
。因此,对于这个初始条件,[0,0]
是所有时间的正确解。
如果将初始条件更改为其他值,您应该会观察到更有趣的结果。