I'm trying to solve a system of differential equations using Python. I've written an algorithm that uses Euler's Method to do so, and I require a time step of 10^-6 s-1, for 100s. That is 10^8 data points, and the computer returns a MemoryError.
The code I have is:
#!/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)
The _new lists are used to help plotting, so as to avoid overloading Matplotlib. Does anyone have any advice on how I can avoid the memory error while still maintaining the required time-step?
PS as part of the project, it is required that I write my own integrator.
I would suggest not trying to keep every iteration of every variable at each time step in memory. You should simply have a 'current' and 'next' version of each variable, update them each time step, and then 'save' the state every 1,000,000 time steps or so. Try something like this:
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
...