如何在python中用矢量化=‘True’实现'solve_ivp‘

Question

我一直试图用solve_ivp来求解一组微分方程。系统的Jacobian矩阵是A，如下所示。我想启用选项vectorized='True'，但不幸的是，我不知道如何修改当前的代码以向量化Jacobian矩阵A。有人知道如何做到这一点吗？

# imports
import numpy as np
import scipy.sparse as sp
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# grid sizing
R=0.05  #sphere radius
N=1000#number of points
D=0.00002 #diffusion coefficient
k=10  # Arrhenius 
Cs=1.0 # Boundary concentration
C0=0.0 # Initial concentration
time_constant=R**2.0/D

dr=R/(N-1)

# Algebra simplification
a=D/dr**2

Init_conc=np.linspace(0,0,N)
B=np.zeros(N)
B[N-1]=Cs*(a+a/(N-1))
#
e1 = np.ones(N)
e2 = np.ones(N)
e3 = np.ones(N)
#
#
#
e1[0]=-k-6*a
e1[1:]=-k-2*a
# 
#
e2[1]=6*a
for i in range(2,N) :
 e2[i]=a+a/(i-1)
#
#
#
for i in range (0,N-1) :
 e3[i]=a-a/(i+1)


A = sp.spdiags([e3,e1,e2],[-1,0,1],N,N,format="csc")


def dc_dt(t,C)   :
    dc=A.dot(C)+B
    return dc

# Solving the system, I want to implement the same thing with vectorized='True'
OutputTimes=np.linspace(0,0.2*time_constant,100)
ans=solve_ivp(dc_dt,(0,0.2*time_constant),Init_conc,method='RK45',t_eval=OutputTimes,vectorized='False')
print (ans)

Answer 1

请看一下this answer，解释是彻底的。有关您的代码，请参见下面更新的代码段和图。并不是很明显，vectorize正在提供任何提速。但是，为关键字jac提供A会产生不同的效果。但我想，只有当A是常数时，它才是有效的？

# imports
import numpy as np
import scipy.sparse as sp
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt  # noqa


def dc_dt(t, C):
    print(C.shape)
    if len(C.shape) == 1:
        return np.squeeze(A.dot(C)) + B
    else:
        return A.dot(C) + np.transpose(np.tile(B, (C.shape[1], 1)))
    # return np.squeeze(A.dot(C)) + B


# grid sizing
R = 0.05  # sphere radius
N = 1000  # number of points
D = 0.00002  # diffusion coefficient
k = 10  # Arrhenius
Cs = 1.0  # Boundary concentration
C0 = 0.0  # Initial concentration
time_constant = R**2.0 / D
dr = R / (N - 1)
# Algebra simplification
a = D / dr**2
Init_conc = np.repeat(0, N)
B = np.zeros(N)
B[-1] = Cs * (a + a / (N - 1))
e1 = np.ones(N)
e2 = np.ones(N)
e3 = np.ones(N)
e1[0] = -k - 6 * a
e1[1:] = -k - 2 * a
e2[1] = 6 * a
for i in range(2, N):
    e2[i] = a + a / (i - 1)
for i in range(0, N - 1):
    e3[i] = a - a / (i + 1)
A = sp.spdiags([e3, e1, e2], [-1, 0, 1], N, N, format="csc")
# Solving the system, I want to implement the same thing with vectorized='True'
OutputTimes = np.linspace(0, 0.2 * time_constant, 10000)
ans = solve_ivp(dc_dt, (0, 0.2 * time_constant), Init_conc,
                method='BDF', t_eval=OutputTimes, jac=A, vectorized=True)
plt.plot(np.arange(N), ans.y[:, 0])
plt.plot(np.arange(N), ans.y[:, 1])
plt.plot(np.arange(N), ans.y[:, 10])
plt.plot(np.arange(N), ans.y[:, 20])
plt.plot(np.arange(N), ans.y[:, 50])
plt.plot(np.arange(N), ans.y[:, -1])
plt.show()

​

​