修复积分器

Question

我正在尝试使用 Python 求解微分方程组。我已经编写了一个使用欧拉方法的算法来执行此操作，并且我需要 10^-6 s-1 的时间步长，持续 100 秒。那是 10^8 个数据点，计算机 returns 一个 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

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

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

Answer 2

我建议不要尝试将每个变量在每个时间步长的每次迭代都保存在内存中。您应该简单地拥有每个变量的 'current' 和 'next' 版本，每个时间步更新它们，然后 'save' 每 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
    ...

修复积分器

Fixing an integrator

python

out-of-memory

integrator

differential-equations