python - 如何在python中使用vectorized='True'实现'solve_ivp'

Question

我一直在尝试使用solve_ivp. 系统的雅可比矩阵是 A，如下所示。我想启用该选项vectorized='True'，但不幸的是我不知道如何修改当前代码以矢量化雅可比矩阵 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)

标签： pythonscipy

Answer 1

请看这个答案，解释很透彻。特别是对于您的代码，请参阅下面的更新片段和图。vectorize提供任何加速并不明显。但是，为关键字提供 Ajac会有所不同。但我想它只有在 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()