CODEWITHSUNDEEP
pythonscipysympy
Tbone Willsone
Tbone Willsone

Reputation: 23

Generating and Solving Simultaneous ODE's in Python

I'm relatively new to Python, and am encountering some issues in writing a piece of code that generates and then solves a system of differential equations.

My approach to doing this was to create a set of variables and coefficients, (x0, x1, ..., xn) and (c0, c1 ,..., cn) repsectively, in a list with the function var(). Then the equations are constructed in EOM1(). The initial conditions, along with the set of equations, are all put together in EOM2() and solved using odeint.

Currently the code below runs, albeit not efficiently the reason for which I believe is because odeint runs through all the code with each interaction (that's something else I need to fix but isn't the main problem!).

import sympy as sy
from scipy.integrate import odeint

n=2
cn0list = [0.01, 0.05]
xn0list = [0.01, 0.01]

def var():
    xnlist=[]
    cnlist=[]
    for i in range(n+1):
        xnlist.append('x{0}'.format(i))
        cnlist.append('c{0}'.format(i))
    return xnlist, cnlist

def EOM1():
    drdtlist=[]
    for i in range(n):
        cn1=sy.Symbol(var()[1][i])
        xn0=sy.Symbol(var()[0][i])
        xn1=sy.Symbol(var()[0][i+1])
        eom=cn1*xn0*(1.0-xn1)-cn1*xn1-xn1
        drdtlist.append(eom)
    xi=sy.Symbol(var()[0][0])
    xf=sy.Symbol(var()[0][n])
    drdtlist[n-1]=drdtlist[n-1].subs(xf,xi)
    return drdtlist

def EOM2(xn, t, cn):
    x0, x1 = xn
    c0, c1 = cn
    f = EOM1()
    output = []
    for part in f:
        output.append(part.evalf(subs={'x0':x0, 'x1':x1, 'c0':c0, 'c1':c1}))
    return output

abserr = 1.0e-6
relerr = 1.0e-4
stoptime = 10.0
numpoints = 20

t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]

wsol = odeint(EOM2, xn0list, t, args=(cn0list,), atol=abserr, rtol=relerr)

My problem is that I had difficulty getting Python to treat the variables generated by Sympy appropriately. I got around this with the line

output.append(part.evalf(subs={'x0':x0, 'x1':x1, 'c0':c0, 'c1':c1}))

in EOM2(). Unfortunately, I do not know how to generalize this line to a list of variables from x0 to xn, and from c0 to cn. The same applies to the earlier line in EOM2(),

    x0, x1 = xn
    c0, c1 = cn

In other words I set n to an arbitrary number, is there a way for Python to interpret each element as it does with the ones I manually entered above? I have tried the following

output.append(part.evalf(subs={'x{0}'.format(j):var(n)[0][j], 'c{0}'.format(j):var(n)[1][j]}))

yet this yields the error that led me to use evalf in the first place,

TypeError: can't convert expression to float

Is there any way do what I want to, generate a set of n equations which are then solved with odeint?

Upvotes: 2

Views: 416

Answers (2)

Bjoern Dahlgren
Bjoern Dahlgren

Reputation: 930

Instead of using evalf you want to look into using sympy.lambdify to generate a callback for use with SciPy. You will need to create a function with the expected signature of odeint, e.g.:

y, params = sym.symbols('y:3'), sym.symbols('kf kb')
ydot = rhs(y, p=params)
f = sym.lambdify((y, t) + params, ydot)
yout = odeint(f, y0, tout, param_values)

We gave a tutorial on (among other things) how to use lambdify with odeint at the SciPy 2017 conference, the material is available here: http://www.sympy.org/scipy-2017-codegen-tutorial/

If you are open to use an external library to handle the function signatures of external solvers you may be interested in a library I've authored: pyodesys

Upvotes: 3

user6655984
user6655984

Reputation:

If I understand correctly, you want to make an arbitrary number of substitutions in a SymPy expression. This is how it can be done: 

n = 10
syms = sy.symbols('x0:{}'.format(n))   # an array of n symbols
expr = sum(syms)                       # some expression with those symbols
floats = [1/(j+1) for j in range(n)]   # numbers to put in 
expr.subs({symbol: value for symbol, value in zip(syms, floats)})

The result of subs is a float in this case (no evalf needed).

Note that the function symbols can directly create any number of symbols for you, via the colon notation. No need for a loop.

Upvotes: 1

Related Questions