如何解决记忆误差并加快Sympy中的方法？

Question

我正在按照稳定性，准确性和其他数值方法来控制runge -kutta方法，以防止过渡超越切换点。

当步骤很小（通过准确性或其他方法）时，整合持续时间很长，并且我遇到了内存错误，

line 712, in _getitem_RepMatrix
    return self._rep.getitem_sympy(index_(i), index_(j))
TypeError: 'slice' object cannot be interpreted as an integer

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
line 157, in extract
    row2[j2] = row1_j1
MemoryError

因此，当我需要大约1E-5或1E-6及以上的准确性时，会发生这种情况。

然后，方法是由方法构建的

def rkf2stepcontrol(
        matrix_size,
        dydx, A, x, b, dxdt, hn, t, tout, eps, predicate_func, MatrixForYacobian, hstabilitygetting, halgebraic):

        k1 = zeros(matrix_size, 1)
        k2 = zeros(matrix_size, 1)
        k3 = zeros(matrix_size, 1)
        k2k1norm = zeros(matrix_size, 1)
        while ((predicate_func == False) or (tout - t > eps)):

            print('Current time', t)
            b = t
            iterations += 1
            print("I am in point ", x)

            xprint_s = xprint_s.col_insert(1, Matrix([x]))

            k1 = hn * dxdt(A, x, b)

            print(k1, 'k1')


            for i in (range(matrix_size)):
                k2[i, 0] = hn * dxdt(A, x + k1[i, 0] * sympy.ones(*x.shape), b)[i, 0]


            xold = x


            for i in (range(matrix_size)):
                x[i, 0] = x[i, 0] + 0.5 * (k1[i, 0] + k2[i, 0])

，通过精确求解方程式0.5” K1 -K2”≤EPS来控制步骤，其中K2和K1-方法步骤和EPS所需的精度。并通过获得最大值（ABS（Jacobian Matrix.EigenValues）来控制稳定性，

所有值将显示在图上，因此我必须包含它们。

Answer 1

切换到numpy会速度提高，但如果您使用大约1E-5或1E-6的精度，则仍然太长了。使用1E-4，整合2.37秒需要55分钟，因此大约一个小时才能制造2.5秒。而且，您只需要按照0.5“ K1-K2”≤EPS来重新计算精度。

代码重写为
...

k1 = numpy.zeros((matrix_size, 1))
k2 = numpy.zeros((matrix_size, 1))
k3 = numpy.zeros((matrix_size, 1))
k2k1norm = numpy.zeros((matrix_size, 1))
xprint_s = numpy.c_[xprint_s, x]
k1 = hn * dxdt(A, x, b)
for i in (range(matrix_size)):
        k2[i, 0] = hn * dxdt(A, x + k1[i, 0] * numpy.ones(x.shape,dtype=float,order='C'), b)[i, 0]



for j in range(matrix_size):            
    k2k1norm[j, 0] = float(dxdt(A, x + k1[i, 0] * numpy.ones(x.shape,dtype=float,order='C'), b)[j, 0]) - \
                                 float(dxdt(A, x, b)[j, 0])
...
return iterations, endtime, xprint_s

并制作情节

yprint_s_rkeiler_angle_step2 = numpy.delete(yprint_s_rkeiler_angle_step2, 0, -1) 
argument1 = list((yprint_s_rkeiler_angle_step2[0]))
argument2 = list((yprint_s_rkeiler_angle_step2[1]))
title('Runge Kutta Esposito Based on Eiler')
plot(t, argument1, '-o', linewidth=2)
plot(t, argument2, '-o', linewidth=2)
legend(["y", "Vy"], loc ="upper right")
print()
ylabel("argument")
xlabel("t")
grid(True)
show() # display