Исправление интегратора - PullRequest
Новая
       19

Исправление интегратора

0 голосов
/ 05 марта 2019

Я пытаюсь решить систему дифференциальных уравнений, используя Python. Я написал алгоритм, который использует метод Эйлера для этого, и мне требуется шаг времени 10 ^ -6 с-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 В рамках проекта требуется написать собственный интегратор.

Ответы [ 2 ]

1 голос
/ 05 марта 2019

Я бы предложил не пытаться хранить каждую итерацию каждой переменной на каждом временном шаге в памяти.Вы должны просто иметь «текущую» и «следующую» версию каждой переменной, обновлять их каждый раз, а затем «сохранять» состояние каждые 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
    ...
0 голосов
/ 05 марта 2019

Вы создаете 8 массивов размером 10 ^ 8. Если каждая запись в массиве является байтом, то это составляет 100 МБ на массив. Даже если вы попытаетесь записать это в файл, это будет довольно большой файл.

Я бы посоветовал вам уменьшить количество точек данных.

...