Python：solve_ivp 的内存分配问题

Question

我在使用 ODE 求解器 solve_ivp 时遇到似乎与内存分配有关的问题。

虽然缩放比例（带有网格点的执行时间）似乎给出了预期的线性结果（一旦我使用了雅可比矩阵的稀疏性），但当我的数组元素数超过 1000 万时，一切都会崩溃。它不应该那样，因为原则上它们应该需要少于 0.5GB 的 RAM 而我有 32GB ...

这里是一个最小的例子，它解决了简单的方程式

y'=1

，但又重现了这个问题。特别是，它运行良好，达到 10^7 点，再多一点就坏了。在我的系统上执行需要 ~O(100) 秒。

import os
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import scipy
from scipy.integrate import solve_ivp
import time
from time import perf_counter
from scipy.sparse import dia_array


def time_integration( Num ):

    t_grid = np.linspace( start=0., stop=1., num=Num )
    data = [np.ones( Num, dtype=np.float32 )]
    jac_sparsity = dia_array( (data, [0]), shape=(Num, Num), dtype=np.float64 )

    def fun_rhs( t, y ):
        dfdt = np.ones( len(t_grid) )
        return dfdt

    start_time = perf_counter()
    output = scipy.integrate.solve_ivp( fun_rhs, [ min(t_grid), max(t_grid) ], np.zeros( len(t_grid) ), 
                                   method='BDF', 
                                   t_eval=None, 
                                   jac_sparsity=jac_sparsity, 
                                   dense_output=False, 
                                   events=None,
                                   vectorized=False,
                                   args=None, atol=1.e-2, rtol=1.e-2 )
    stop_time = perf_counter()
    return stop_time - start_time


start_log = 1.
stop_log = 7.2   ## !!!!! if 6. (or even 7.) it works fine !!!!!
timeSteps_array = np.logspace(start=start_log, stop=stop_log, num=int( stop_log-start_log ) + 1, dtype=int)
print(f'array of time steps (len = {len(timeSteps_array)}): {timeSteps_array}')
print('')


start_time_script = perf_counter()
runtime_array = np.zeros( len(timeSteps_array) )
runtime_array = [time_integration( it ) for it in timeSteps_array]
stop_time_script = perf_counter()


print(f'runtime for each Nt = {runtime_array} sec')
print(f'the script took {stop_time_script - start_time_script} seconds to run.')
plt.loglog(timeSteps_array, runtime_array, lw=2., color='blue', label='num test')
plt.loglog(timeSteps_array, time_integration( timeSteps_array[1] ) * timeSteps_array/timeSteps_array[1], lw=1.5, ls='--', color='red', label='linear')
plt.xlabel('Total number of points', fontsize=13)
plt.ylabel('Time [sec]', fontsize=13)
plt.legend(frameon=False, fontsize=12)
plt.savefig('TimeScaling.pdf',format='pdf',bbox_inches='tight', dpi=200)
plt.show()

补充一点应该有用，我得到的错误是：

RuntimeError: SUPERLU_MALLOC fails for buf in intCalloc() at line 159 in file /private/tmp/scipy-20230221-70667-mabqq4/scipy-1.10.1/scipy/sparse/linalg/_dsolve/SuperLU/SRC/memory.c