CODEWITHSUNDEEP
pythonscipy
Laura Holewa
Laura Holewa

Reputation: 121

How can I modify this differential equation solver to solve for a large number of variables?

I would like to modify the following code so that it can solve for thousands of variables with thousands of couple differential equations. The problem is that I would like to be able to import the variable list (y), ideally as a numpy array, but acceptably as a regular list. The list will be massive, so I don't want to define it in the function. If I do something like define a='n','c1',... and then set a=y, python complains that the variables are the wrong type. If I define a= n, c1,....python complains that I haven't defined the variables. If you have any advice on how to import the very large variable list for this function as a list or numpy array that would be great. 

import numpy as np
from scipy.integrate import odeint


def kinetics(y,t,b1,b2):
    n,c1=y
    dydt=[(b1*n)+(b2*c1),(b1*n)-(c1)]
    return dydt

b1=0.00662888
b2=0.000239997

n0=1
c1_0=1

y0=[n0,c1_0]
t = np.linspace(0, 10, 10)
sol = odeint(kinetics, y0, t, args=(b1,b2))

print(sol)

Upvotes: 0

Views: 257

Answers (1)

Matthew Flamm
Matthew Flamm

Reputation: 620

This sort of renaming and use of individual python objects for each step is likely to be slow, especially if you can use numpy array vectorization for parts of the computation. Note that your current code can be reformulated into a linear algebra problem which could be solved very efficiently using numpy, assuming this is true for the rest of the ODE structure. I warn against the path you want to take.

But, assuming that nothing can be vectorized in this way, and you want to continue down this path, you could define a class to store all your data and create an instance of said class with the current state of y at each function call. You could store the class definition in a separate python module to keep the problem definition tidy if needed.

If using python 3:

import numpy as np
from scipy.integrate import odeint


class vars:
    def __init__(self, *args):
        (self.n, self.c1) = args

def kinetics(y,t,b1,b2):
    v = vars(*y)
    dydt=[
        (b1 * v.n) + (b2 * v.c1),
        (b1 * v.n) - (v.c1)
    ]
    return dydt

b1=0.00662888
b2=0.000239997

n0=1
c1_0=1

y0=[n0,c1_0]
t = np.linspace(0, 10, 10)
sol = odeint(kinetics, y0, t, args=(b1,b2))

print(sol)

Upvotes: 1

Related Questions