我能想到的一些事情：

• 使用稀疏矩阵来加速Matrix1和的计算Matrix2dZdY
• 使用更大的积分容差reltol和abstol，特别是如果您正在寻找稳态解决方案和/或不需要精确解析系统瞬态动态。
• 您正在使用 ode15s 求解，它是一个隐式积分器，您应该自己提供 ODE 函数的雅可比行列式，因为在这里看起来很简单。这是通过关键字“Jacobian”完成的。
• 如果 Jacobian 太重而无法导出（可能是更复杂的问题），您至少可以通过关键字“JPattern”指定其稀疏模式。如果雅可比是块对角线，则求解器将能够通过有限差分计算它，而调用次数比考虑密集雅可比（否则为默认值）要少得多。
• 最后，对函数进行矢量化并将关键字“矢量化”设置为 True。这样，求解器可以一次调用具有多个输入的 ODE 函数并减少一些开销。这对于之前的雅可比计算特别有用。

在这里，我认为代码可以很容易地通过简单地将等式重写为（+ 或 - 对元素操作的一些调整）进行矢量化：

dY(1,:) = ...
dY(end,:) = ...
dY(2:end-1,:) = ... 

哦，当然，如果您的问题是非刚性 ODE，只需使用像 之类的显式求解器ode45，它们会快得多。除非您想要达到稳态并且初始瞬态很长，在这种情况下，具有较大误差容限的隐式方法可能能够更快地达到稳态，而无需精细解决瞬态。

编辑：如果我误解了您的答案并且您实际上想要运行 2400 个不同的模拟，那么将它们作为单个“大”模拟运行可能不是一个好主意，因为时间步长将受到“最难”解决方案的限制。在这种情况下，另一种选择是使用并行 for 循环（使用parfor）并在多个内核上分别解决每种情况。

编辑2021 年 2 月 15 日关于稀疏 aptternodenumjac的形成，您可以这样做：

f0 = ode(0,y0); % value of your ODE function
J = odenumjac(ode,{0 y0}, f0); ! Computes the Jacobian, assuming it is dense
sparsity_pattern = sparse(Jac~=0.);

然后，您可以在调用 ODE 求解器时将其指定'JPattern', sparsity_pattern为附加选项。Matlab 将利用这种模式的知识通过有限差分更快地逼近雅可比行列式。请注意，在您的情况下，您的雅可比很容易推导，因此您最好通过选项通过函数提供雅可比，Jacobian或者，如果它是常数，则作为常数雅可比 via JConstant。

关于odenumjac近似雅可比的方式，最好的办法是在互联网上查找“雅可比有限差分”以找到更深入的解释，但这是主要原理。你 ODE 函数$f$接受一个输入（状态向量$x$) 大小$N$. 一些 ODE 求解器（通常是隐式求解器）需要 Jacobian NxN$ 矩阵）的知识才能计算解。$\partial_x f` ($

行是的第个分量相对于的每个分量的偏导数的向量。相反，列是的所有分量相对于个分量的偏导数的向量。该列可以通过以下一阶有限差分公式近似：$i$$\partial_x f$$i$$f$$x$$i$$\partial_x f$$f$$i$$y$

上进行循环，可以近似整个雅可比行列式。这将需要次评估：一次不受干扰地计算，而次评估来计算。$i$$N+1$$f$$f(x)$$N$$f(x+h e_i), i=1..N$

对于稀疏雅可比结构的一般问题（例如离散 PDE），当的分量没有“重叠”影响时，稀疏模式的知识可用于一次执行多个扰动。您可以查看文章“关于稀疏雅可比矩阵的估计”以获取更多信息：http ://www.ii.uib.no/forskningsgrupper/opt/forskning/papers/extern/cpr.pdf$x$$f$

发布的简单示例并不僵硬，但无论如何我还是将它放在 Julia 中进行展示。我将它修改为一个包含 2 个 PDE 的系统，N=1200并进入 10 毫秒

using ModelingToolkit, LinearAlgebra, BenchmarkTools, DifferentialEquations

# Setup matrices
N = 1200
mat1 = hcat(zeros(N),Tridiagonal(ones(N-1),-2*ones(N),ones(N-1)),zeros(N))
mat1[1,1] = 1
mat1[end,end] = 1

mat2 = hcat(zeros(N),Bidiagonal(-1*ones(N),ones(N-1), :L),zeros(N))
mat2[1,1] = 1

# Define the ODE

function f(du,u,p,t)
    y = @view u[:,1]
    z = @view u[:,2]
    dy = @view du[:,1]
    dz = @view du[:,2]

    dy[1] = 0
    dy[2:end-1] .= mat1*y + mat2*y
    dy[end] = 2*(y[end-1] - y[end])

    dz[1] = 0
    dz[2:end-1] .= mat1*z + mat2*z
    dz[end] = 2*(z[end-1] - z[end])
end
u0 = zeros(N+2,2); u0[1,:] .= 1
tspan = (0.0,5.0)
prob = ODEProblem(f, u0, tspan)

# Accelerate it
sys = modelingtoolkitize(prob)
fastprob = ODEProblem(sys, u0, tspan, jac=true, sparse=true)

# Non-Stiff Choice
# 12.908 ms (1073 allocations: 9.58 MiB)
@btime solve(fastprob,Tsit5(),saveat=0.01,abstol=1e-10,reltol=1e-9)

# Unaccelerated Non-Stiff Choice
# 21.300 ms (9658 allocations: 36.36 MiB)
@btime solve(prob,Tsit5(),saveat=0.01,abstol=1e-10,reltol=1e-9)

# Stiff Choice
# 479.284 ms (21441 allocations: 619.75 MiB)
@btime solve(fastprob,TRBDF2(),saveat=0.01,abstol=1e-10,reltol=1e-9)

# Unaccelerated Stiff Choice
# 15.849 s (1381204 allocations: 385.17 MiB)
@btime solve(prob,TRBDF2(),saveat=0.01,abstol=1e-10,reltol=1e-9)

您可能可以将其修改为原始问题。