有没有一种方法可以轻松地在整个点网格上集成一组微分方程？
Is there a way to easily integrate a set of differential equations over a full grid of points?
问题是我希望能够一次对从网格的每个点开始的微分方程进行积分，而不必为每个坐标循环遍历 scipy 积分器。 （我确定有一个简单的方法）
作为代码的背景，我试图解决每个特定周期交替变化速度方向的库埃特通量的轨迹，这是一个众所周知的产生混沌的动力系统。作为与 scipy 集成的一部分以及我对 numpy 的 meshgrid 函数的使用，我认为其余代码并不重要。
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, writers
from scipy.integrate import solve_ivp
start_T = 100
L = 1
V = 1
total_run_time = 10*3
grid_points = 10
T_list = np.arange(start_T, 1, -1)
x = np.linspace(0, L, grid_points)
y = np.linspace(0, L, grid_points)
X, Y = np.meshgrid(x, y)
condition = True
totals = np.zeros((start_T, total_run_time, 2))
alphas = np.zeros(start_T)
i = 0
for T in T_list:
alphas[i] = L / (V * T)
solution = np.array([X, Y])
for steps in range(int(total_run_time/T)):
t = steps*T
if condition:
def eq(t, x):
return V * np.sin(2 * np.pi * x[1] / L), 0.0
condition = False
else:
def eq(t, x):
return 0.0, V * np.sin(2 * np.pi * x[1] / L)
condition = True
time_steps = np.arange(t, t + T)
xt = solve_ivp(eq, time_steps, solution)
solution = np.array([xt.y[0], xt.y[1]])
totals[i][t: t + T][0] = solution[0]
totals[i][t: t + T][1] = solution[1]
i += 1
np.save('alphas.npy', alphas)
np.save('totals.npy', totals)
给出的错误是：
ValueError: y0 must be 1-dimensional.
而它来自scipy的'solve_ivp'函数，因为它不接受numpy函数meshgrid的格式。我知道我可以 运行 一些循环并克服它，但我假设必须有一个 'good' 方法来使用 numpy 和 scipy 来完成它。我也接受对其余代码的建议。
TL;DR
我不认为你可以“一次对从网格的每个点开始的微分方程进行积分”。
MWE
请尝试提供 MWE 来重现您的问题，就像您说的：“我认为其余代码并不重要”，这会让人们更难理解您的问题.
了解如何与求解器对话
在回答你的问题之前，有几件事似乎被误解了：
- 通过定义
time_steps = np.arange(t, t + T) 然后调用 solve_ivp(eq, time_steps, solution) ：solve_ivp 的第二个参数是您想要解决方案的时间跨度，即“开始”和“停止”作为 2-uple 的时间。这里你的 time_steps 是 30 长（对于第一个循环），所以我可能会用 (t, t+T) 替换它。 Look for t_span in the doc.
- 据我了解，您似乎想要控制数值分辨率的每次迭代：这不是
solve_ivp 的工作方式。此外，我认为您想在每次迭代时切换函数“eq”。由于您必须传递等式的“右侧”，因此您需要将此行为 包装在 一个函数中。它不会工作（见后面）但就概念而言是这样的：
def RHS(t, x):
# unwrap your variables, condition is like an additional variable of your problem,
# with a very simple differential equation
x0, x1, condition = x
# compute new results for x0 and x1
if condition:
x0_out, x1_out = V * np.sin(2 * np.pi * x[1] / L), 0.0
else:
x0_out, x1_out = 0.0, V * np.sin(2 * np.pi * x[1] / L)
# compute new result for condition
condition_out = not(condition)
return [x0_out, x1_out, condition_out]
这是行不通的，因为条件的演化不满足derivation/continuity的一些数学性质。所以 condition 就像一个参数化模型的布尔开关，我们可以使用 global 来控制这个布尔值的状态：
condition = True
def RHS_eq(t, y):
global condition
x0, x1 = y
# compute new results for x0 and x1
if condition:
x0_out, x1_out = V * np.sin(2 * np.pi * x1 / L), 0.0
else:
x0_out, x1_out = 0.0, V * np.sin(2 * np.pi * x1 / L)
# update condition
condition = 0 if condition==1 else 1
return [x0_out, x1_out]
- 最后，这是您在 post 中提到的
ValueError：您定义 solution = np.array([X, Y])，它实际上是初始条件，应该是“y0: array_like, shape (n,)”，其中n 是问题的变量数（在 [x0_out, x1_out] 的情况下为 2）
单个初始条件的 MWE
综上所述，让我们从单个起点 (0.5,0.5) 的简单 MWE 开始，这样我们就可以清楚地了解如何使用求解器：
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# initial conditions for x0, x1, and condition
initial = [0.5, 0.5]
condition = True
# time span
t_span = (0, 100)
# constants
V = 1
L = 1
# define the "model", ie the set of equations of t
def RHS_eq(t, y):
global condition
x0, x1 = y
# compute new results for x0 and x1
if condition:
x0_out, x1_out = V * np.sin(2 * np.pi * x1 / L), 0.0
else:
x0_out, x1_out = 0.0, V * np.sin(2 * np.pi * x1 / L)
# update condition
condition = 0 if condition==1 else 1
return [x0_out, x1_out]
solution = solve_ivp(RHS_eq, # Right Hand Side of the equation(s)
t_span, # time span, a 2-uple
initial, # initial conditions
)
fig, ax = plt.subplots()
ax.plot(solution.t,
solution.y[0],
label="x0")
ax.plot(solution.t,
solution.y[1],
label="x1")
ax.legend()
最终答案
现在，我们想要做的是完全相同的事情，但对于不同的初始条件，据我所知，我们不能：再次引用文档
y0 : array_like, shape (n,) : Initial state. 。求解器的初始条件只允许一个起点向量。
所以回答最初的问题：我不认为你可以“立即对从网格的每个点开始的微分方程进行积分”。
是的，您可以通过多种方式做到这一点。如果建议这样做，问题仍然存在。
要实现一个普遍可用的 ODE 积分器，它需要从模型中抽象出来。大多数实现通过让状态 space 成为一个平面数组向量 space 来做到这一点，有些实现允许将向量 space 引擎作为参数传递，因此结构化向量 spaces可以使用。 scipy 集成商不是这种类型。
因此您需要将状态转换为积分器的平面向量，然后再转换回模型的结构化状态。
def encode(X,Y): return np.concatenate([X.flatten(),Y.flatten()])
def decode(U): return U.reshape([2,grid_points,grid_points])
然后你可以将ODE函数实现为
def eq(t,U):
X,Y = decode(U)
Vec = V * np.sin(2 * np.pi * x[1] / L)
if int(t/T)%2==0:
return encode(Vec, np.zeros(Vec.shape))
else:
return encode(np.zeros(Vec.shape), Vec)
初始值
U0 = encode(X,Y)
那么可以直接在整个时间跨度上进行整合。
为什么这可能不是一个好主意：分别考虑每个网格点及其轨迹，对于给定的误差水平，每个轨迹都有自己的自适应时间步长序列。在同时整合所有轨迹时，调整后的步长是给定时间所有轨迹中的最小步长。因此，虽然单个轨迹在时间步长稀疏的长间隔中可能只有很短的间隔和非常小的步长，但这些轨迹可能会在整体中重叠，从而导致到处都是非常小的步长。
如果您超出了测试阶段，请切换到编译程度更高的求解器实现，odeint 是带有包装器的 Fortran 代码，所以是半个解决方案。 JITcode 转换为 C 代码并链接到 odeint 后面的已编译求解器。离开 python 你会得到日晷，julia-lang 的 diffeq 模块，或者 boost::odeint.
问题是我希望能够一次对从网格的每个点开始的微分方程进行积分，而不必为每个坐标循环遍历 scipy 积分器。 （我确定有一个简单的方法）
作为代码的背景，我试图解决每个特定周期交替变化速度方向的库埃特通量的轨迹，这是一个众所周知的产生混沌的动力系统。作为与 scipy 集成的一部分以及我对 numpy 的 meshgrid 函数的使用，我认为其余代码并不重要。
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, writers
from scipy.integrate import solve_ivp
start_T = 100
L = 1
V = 1
total_run_time = 10*3
grid_points = 10
T_list = np.arange(start_T, 1, -1)
x = np.linspace(0, L, grid_points)
y = np.linspace(0, L, grid_points)
X, Y = np.meshgrid(x, y)
condition = True
totals = np.zeros((start_T, total_run_time, 2))
alphas = np.zeros(start_T)
i = 0
for T in T_list:
alphas[i] = L / (V * T)
solution = np.array([X, Y])
for steps in range(int(total_run_time/T)):
t = steps*T
if condition:
def eq(t, x):
return V * np.sin(2 * np.pi * x[1] / L), 0.0
condition = False
else:
def eq(t, x):
return 0.0, V * np.sin(2 * np.pi * x[1] / L)
condition = True
time_steps = np.arange(t, t + T)
xt = solve_ivp(eq, time_steps, solution)
solution = np.array([xt.y[0], xt.y[1]])
totals[i][t: t + T][0] = solution[0]
totals[i][t: t + T][1] = solution[1]
i += 1
np.save('alphas.npy', alphas)
np.save('totals.npy', totals)
给出的错误是：
ValueError: y0 must be 1-dimensional.
而它来自scipy的'solve_ivp'函数，因为它不接受numpy函数meshgrid的格式。我知道我可以 运行 一些循环并克服它，但我假设必须有一个 'good' 方法来使用 numpy 和 scipy 来完成它。我也接受对其余代码的建议。
TL;DR
我不认为你可以“一次对从网格的每个点开始的微分方程进行积分”。
MWE
请尝试提供 MWE 来重现您的问题，就像您说的：“我认为其余代码并不重要”，这会让人们更难理解您的问题.
了解如何与求解器对话
在回答你的问题之前，有几件事似乎被误解了：
- 通过定义
time_steps = np.arange(t, t + T)然后调用solve_ivp(eq, time_steps, solution)：solve_ivp的第二个参数是您想要解决方案的时间跨度，即“开始”和“停止”作为 2-uple 的时间。这里你的time_steps是 30 长（对于第一个循环），所以我可能会用(t, t+T)替换它。 Look fort_spanin the doc. - 据我了解，您似乎想要控制数值分辨率的每次迭代：这不是
solve_ivp的工作方式。此外，我认为您想在每次迭代时切换函数“eq”。由于您必须传递等式的“右侧”，因此您需要将此行为 包装在 一个函数中。它不会工作（见后面）但就概念而言是这样的：
def RHS(t, x):
# unwrap your variables, condition is like an additional variable of your problem,
# with a very simple differential equation
x0, x1, condition = x
# compute new results for x0 and x1
if condition:
x0_out, x1_out = V * np.sin(2 * np.pi * x[1] / L), 0.0
else:
x0_out, x1_out = 0.0, V * np.sin(2 * np.pi * x[1] / L)
# compute new result for condition
condition_out = not(condition)
return [x0_out, x1_out, condition_out]
这是行不通的，因为条件的演化不满足derivation/continuity的一些数学性质。所以 condition 就像一个参数化模型的布尔开关，我们可以使用 global 来控制这个布尔值的状态：
condition = True
def RHS_eq(t, y):
global condition
x0, x1 = y
# compute new results for x0 and x1
if condition:
x0_out, x1_out = V * np.sin(2 * np.pi * x1 / L), 0.0
else:
x0_out, x1_out = 0.0, V * np.sin(2 * np.pi * x1 / L)
# update condition
condition = 0 if condition==1 else 1
return [x0_out, x1_out]
- 最后，这是您在 post 中提到的
ValueError：您定义solution = np.array([X, Y])，它实际上是初始条件，应该是“y0: array_like, shape (n,)”，其中n 是问题的变量数（在[x0_out, x1_out]的情况下为 2）
单个初始条件的 MWE
综上所述，让我们从单个起点 (0.5,0.5) 的简单 MWE 开始，这样我们就可以清楚地了解如何使用求解器：
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# initial conditions for x0, x1, and condition
initial = [0.5, 0.5]
condition = True
# time span
t_span = (0, 100)
# constants
V = 1
L = 1
# define the "model", ie the set of equations of t
def RHS_eq(t, y):
global condition
x0, x1 = y
# compute new results for x0 and x1
if condition:
x0_out, x1_out = V * np.sin(2 * np.pi * x1 / L), 0.0
else:
x0_out, x1_out = 0.0, V * np.sin(2 * np.pi * x1 / L)
# update condition
condition = 0 if condition==1 else 1
return [x0_out, x1_out]
solution = solve_ivp(RHS_eq, # Right Hand Side of the equation(s)
t_span, # time span, a 2-uple
initial, # initial conditions
)
fig, ax = plt.subplots()
ax.plot(solution.t,
solution.y[0],
label="x0")
ax.plot(solution.t,
solution.y[1],
label="x1")
ax.legend()
最终答案
现在，我们想要做的是完全相同的事情，但对于不同的初始条件，据我所知，我们不能：再次引用文档
y0 : array_like, shape (n,) : Initial state. 。求解器的初始条件只允许一个起点向量。
所以回答最初的问题：我不认为你可以“立即对从网格的每个点开始的微分方程进行积分”。
是的，您可以通过多种方式做到这一点。如果建议这样做，问题仍然存在。
要实现一个普遍可用的 ODE 积分器，它需要从模型中抽象出来。大多数实现通过让状态 space 成为一个平面数组向量 space 来做到这一点，有些实现允许将向量 space 引擎作为参数传递，因此结构化向量 spaces可以使用。 scipy 集成商不是这种类型。
因此您需要将状态转换为积分器的平面向量，然后再转换回模型的结构化状态。
def encode(X,Y): return np.concatenate([X.flatten(),Y.flatten()])
def decode(U): return U.reshape([2,grid_points,grid_points])
然后你可以将ODE函数实现为
def eq(t,U):
X,Y = decode(U)
Vec = V * np.sin(2 * np.pi * x[1] / L)
if int(t/T)%2==0:
return encode(Vec, np.zeros(Vec.shape))
else:
return encode(np.zeros(Vec.shape), Vec)
初始值
U0 = encode(X,Y)
那么可以直接在整个时间跨度上进行整合。
为什么这可能不是一个好主意：分别考虑每个网格点及其轨迹，对于给定的误差水平，每个轨迹都有自己的自适应时间步长序列。在同时整合所有轨迹时，调整后的步长是给定时间所有轨迹中的最小步长。因此，虽然单个轨迹在时间步长稀疏的长间隔中可能只有很短的间隔和非常小的步长，但这些轨迹可能会在整体中重叠，从而导致到处都是非常小的步长。
如果您超出了测试阶段，请切换到编译程度更高的求解器实现，odeint 是带有包装器的 Fortran 代码，所以是半个解决方案。 JITcode 转换为 C 代码并链接到 odeint 后面的已编译求解器。离开 python 你会得到日晷，julia-lang 的 diffeq 模块，或者 boost::odeint.