修复积分器答案

【问题标题】：Fixing an integrator修复积分器
【发布时间】：2019-03-05 14:24:07
【问题描述】：

我正在尝试使用 Python 求解微分方程组。我编写了一个使用欧拉方法的算法，我需要 10^-6 s-1 的时间步长，持续 100 秒。那是 10^8 个数据点，计算机返回 MemoryError。

我的代码是：

#!/usr/bin/env python3

import matplotlib.pyplot as plt
import math
import numpy as np

k1 = 1.34
k2 = 1.6E+9
k3 = 8E+3
k4 = 4E+7
k5 = 1

def f_A(A,Y):
    return -k1*A*Y

def f_B(B,X):
    return -k3*X*B

def f_X(X,Y,A,B):
    return k1*A*Y - k2*X*Y + k3*B*X - k4*X*X 

def f_Y(X,Y,Z,A):
    return -k1*A*Y - k2*X*Y + k5*Z

def f_Z(X,Z,B):
    return -k5*Z + k3*B*X

def f_P(X,Y,A):
    return k1*A*Y + k2*X*Y

def f_Q(X):
    return k4*X*X

def Euler(fA,fB,fX,fY,fZ,fP,fQ,t0,tt,n):
    h = (tt - t0) / float(n)

    t = [0]*(n)
    X = [0]*(n)
    Y = [0]*(n)
    Z = [0]*(n)
    P = [0]*(n)
    Q = [0]*(n)
    A = [0]*(n)
    B = [0]*(n)

    t[0] = t0
    X[0] = 10**-9.8
    Y[0] = 10**-6.52
    Z[0] = 10**-7.32
    A[0] = 0.06
    B[0] = 0.06
    P[0] = 0
    Q[0] = 0

    for i in range(1,n):

        t[i] = t0 + i*h

        X[i] = X[i-1] + h*fX(X[i-1],Y[i-1],A[i-1],B[i-1])

        Y[i] = Y[i-1] + h*fY(X[i-1],Y[i-1], Z[i-1], A[i-1])

        Z[i] = Z[i-1] + h*fZ(X[i-1],Z[i-1],B[i-1])

        A[i] = A[i-1] + h*fA(A[i-1],Y[i-1])

        B[i] = B[i-1] + h*fB(B[i-1],X[i-1])

        P[i] = P[i-1] + h*fP(X[i-1],Y[i-1],A[i-1])

        Q[i] = Q[i-1] + h*fQ(X[i-1])

    t_new = t[0::100]
    X_new = X[0::100]
    Y_new = Y[0::100]
    Z_new = Z[0::100]


    plt.figure(figsize=(10, 4))
    plt.yscale('log')
    plt.plot(t_new, X_new, label = 'X')
    plt.plot(t_new, Y_new, label = 'Y')
    plt.plot(t_new, Z_new, label = 'Z')
    plt.xlabel('time / s')
    plt.ylabel('concentration')
    plt.legend()
    plt.show()

t_0 = 0
t_t = 100 
m = 10**8

Euler(f_A,f_B,f_X,f_Y,f_Z,f_P,f_Q,t_0,t_t,m)

_new 列表用于帮助绘图，以避免 Matplotlib 重载。有人对我如何在保持所需时间步长的同时避免内存错误有任何建议吗？

PS 作为项目的一部分，需要我自己编写积分器。

【问题讨论】：

你真的需要绘制 10^8 个数据点吗？尝试为具有 10^8 个元素的列表分配内存可能是导致错误的原因。
不作图测试，为了可视化做一些平均
除非您这样做是为了练习，或者您真的知道自己在做什么，否则请不要自己使用欧拉积分器或程序积分器。那个轮子has already been invented.
这确实是为了运动而伤心...

标签： python out-of-memory differential-equations integrator

【解决方案1】：

我建议不要尝试将每个变量在每个时间步的每次迭代都保留在内存中。您应该简单地拥有每个变量的“当前”和“下一个”版本，在每个时间步更新它们，然后每 1,000,000 个时间步左右“保存”状态。试试这样的：

def Euler(fA,fB,fX,fY,fZ,fP,fQ,t0,tt,n):
    num_samples = 100

    h = (tt - t0) / float(n)

    # initialise variables
    t = t0
    X = 10**-9.8
    ...

    # initialise _samples lists
    t_samples = []
    X_samples = []
    Y_samples = []
    Z_samples = []

    for i in range(1,n):
        # save the state once every (n / num_samples) time steps
        if i % (n / num_samples) == 0:
            t_samples.append(t)
            X_samples.append(X)
            Y_samples.append(Y)
            Z_samples.append(Z)

        # compute the next version of each variable
        t_ = t0 + i*h
        X_ = X + h*fX(X, Y, A, B)
        ...

        # update the variables
        t, X, Y, Z, A, B, P, Q = t_, X_, Y_, Z_, A_, B_, P_, Q_

    # plot using _samples lists
    ...

【讨论】：

【解决方案2】：

您正在创建 8 个大小为 10^8 的数组。如果数组中的每个条目都是一个字节，那么这相当于每个数组 100 MB。即使您尝试将其写入文件，它也将是一个非常大的文件。

我建议你必须减少数据点的数量。

【讨论】：