Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

SymPy/SciPy: solving a system of ordinary differential equations with different variables

I am new to SymPy and Python in general, and I am currently working with Python 2.7 and SymPy 0.7.5 with the objective to: a) read a system of differential equations from a text file b) solve the system

I already read this question and this other question, and they are almost what I am looking for, but I have an additional issue: I do not know in advance the form of the system of equations, so I cannot create the corresponding function using def inside the script, as in this example. The whole thing has to be managed at run-time.

So, here are some snippets of my code. Suppose I have a text file system.txt containing the following:

dx/dt = 0.0387*x - 0.0005*x*y
dy/dt = 0.0036*x*y - 0.1898*y

What I do is:

# imports
import sympy
import scipy
import re as regex

# define all symbols I am going to use
x = sympy.Symbol('x')
y = sympy.Symbol('y')
t = sympy.Symbol('t')

# read the file
systemOfEquations = []
with open("system.txt", "r") as fp :
   for line in fp :
            pattern = regex.compile(r'.+?\s+=\s+(.+?)$')
            expressionString = regex.search(pattern, line) # first match ends in group(1)   
            systemOfEquations.append( sympy.sympify( expressionString.group(1) ) )

At this point, I am stuck with the two symbolic expressions inside the systemOfEquation list. Provided that I can read the initial conditions for the ODE system from another file, in order to use scipy.integrate.odeint, I would have to convert the system into a Python-readable function, something like:

def dX_dt(X, t=0):
return array([ 0.0387*X[0] - 0.0005*X[0]*X[1] ,
              -0.1898*X[1] + 0.0036*X[0]*X[1] ])

Is there a nice way to create this at run-time? For example, write the function to another file and then import the newly created file as a function? (maybe I am being stupid here, but remember that I am relatively new to Python :-D)

I've seen that with sympy.utilities.lambdify.lambdify it's possible to convert a symbolic expression into a lambda function, but I wonder if this can help me...lambdify seems to work with one expression at the time, not with systems.

Thank you in advance for any advice :-)

EDIT:

With minimal modifications, Warren's answer worked flawlessly. I have a list of all symbols inside listOfSymbols; moreover, they appear in the same order as the columns of data X that will be used by odeint. So, the function I used is

def dX_dt(X, t):
    vals = dict()
    for index, s in enumerate(listOfSymbols) :
            if s != time :
                    vals[s] = X[index]
    vals[time] = t
    return [eq.evalf(subs=vals) for eq in systemOfEquations]

I just make an exception for the variable 'time' in my specific problem. Thanks again! :-)

like image 241
Alberto Avatar asked Apr 03 '14 10:04

Alberto


People also ask

How do you solve equations with variables in Python?

To solve the two equations for the two variables x and y , we'll use SymPy's solve() function. The solve() function takes two arguments, a tuple of the equations (eq1, eq2) and a tuple of the variables to solve for (x, y) . The SymPy solution object is a Python dictionary.

How do I isolate a variable in SymPy?

To isolate one variable, call the solve function, and pass that variable as the second argument. The brackets surround a list of all solutions—in this case, just one. That solution is that R = P V T n .

How do you solve differential equations in Python?

Differential equations are solved in Python with the Scipy. integrate package using function odeint or solve_ivp. t: Time points at which the solution should be reported. Additional internal points are often calculated to maintain accuracy of the solution but are not reported.


1 Answers

If you are going to solve the system in the same script that reads the file (so systemOfEquations is available as a global variable), and if the only variables used in systemOfEquations are x, y and possibly t, you could define dX_dt in the same file like this:

def dX_dt(X, t):
    vals = dict(x=X[0], y=X[1], t=t)
    return [eq.evalf(subs=vals) for eq in systemOfEquations]

dX_dt can be used in odeint. In the following ipython session, I have already run the script that creates systemOfEquations and defines dX_dt:

In [31]: odeint(dX_dt, [1,2], np.linspace(0, 1, 5))
Out[31]: 
array([[ 1.        ,  2.        ],
       [ 1.00947534,  1.90904183],
       [ 1.01905178,  1.82223595],
       [ 1.02872997,  1.73939226],
       [ 1.03851059,  1.66032942]]
like image 53
Warren Weckesser Avatar answered Sep 30 '22 18:09

Warren Weckesser