修复集成商

Question

我正在尝试使用Python解决微分方程组。我编写了一个使用Euler方法的算法，我需要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作为项目的一部分，我需要编写自己的集成商。

Answer 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
    ...

Answer 2

您正在创建8个大小为10 ^ 8的数组。如果数组中的每个条目都是一个字节，则每个数组的数量为100 MB。即使你试图将其写入文件，它也会是一个非常大的文件。

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