best way to solve time-dependent series of equations

40 views
Skip to first unread message

John

unread,
May 16, 2015, 3:56:47 PM5/16/15
to sy...@googlegroups.com
Hello:
I am seeking the most elegant way to use Sympy to solve a time-dependent series of equations.  As a simple example, consider a three-compartment model with compartments A, B, and C.  Flow occur from A to B, and from B to C.  The following solves such a system that is not time-dependent:

import sympy as sym
import numpy as np

fractAB = .1
fractBC = .2

A, B, C, flowAB, flowBA, flowBC, flowCB = sym.S('A, B, C, flowAB, flowBA, flowBC, flowCB'.split(',')) 

equations = [
  sym.Eq(A,10),
  sym.Eq(flowAB, A*fractAB),
  sym.Eq(flowBC, B*fractBC),
  sym.Eq(A, flowBA - flowAB),
  sym.Eq(B, flowAB + flowCB - flowBA - flowBC),
  sym.Eq(C, flowBC - flowCB),
  sym.Eq(C, 0)
  ]   
 
res = sym.solve(equations)

My question is, what is the most elegant way to solve the system of equations if fractAB and fractBC are both functions of time (perhaps specified as a numpy array, with each element representing a time step).  Can someone please suggest a snipit of code to do this?  One way would be to hold fractAB and fract BC in a 2xn matrix, M, where n is the number of time steps.  Then the code above could simply be put in a for loop, as in:

for i in range(n):
    fractAB = M[0,i]
    fractBC = M[1,i]
    ... same code as above, with each iteration producing a result from the call to solve.

But there is probably a more elegant way to do this.  Thanks in advance for the suggestions.

John

unread,
May 16, 2015, 7:21:15 PM5/16/15
to sy...@googlegroups.com
I have tried to use dsolve, but when I run my code I get errors. I think I am using the wrong syntax, or setting up the problem wrong.  Consider the even simpler two-compartment model at steady state.  The following code works fine:

import sympy as sym

fractAB = .1
fractBA = .2

A, B, flowAB, storeB, flowBA = sym.S('A, B, storeB, flowAB, flowBA'.split(',')) 

equations = [
  sym.Eq(A,10),
  sym.Eq(flowAB, A*fractAB),
  sym.Eq(A, flowBA - flowAB),
  sym.Eq(B, storeB + flowAB - flowBA),
  sym.Eq(flowBA, B*fractBA)
  ]   
 
res = sym.solve(equations)

To make fractAB and fract BA functions of time, which is what I am after, I tried the following:

import sympy as sym

t, A = sym.S('t, A'.split(','))

B, flowAB, storeB, flowBA, fractAB, fractBA = sym.symbols(
  'B, storeB, flowAB, flowBA, fractAB, fractBA'.split(','), cls=sym.Function) 

equations = (
  sym.Eq(A,10),
  sym.Eq(flowAB(t), A*fractAB(t)),
  sym.Eq(A, flowBA(t) - flowAB(t)),
  sym.Eq(B(t), storeB(t) + flowAB(t) - flowBA(t)),
  sym.Eq(flowBA(t), B(t) * fractBA(t)),
  sym.Eq(fractAB(t), .1*t),
  sym.Eq(fractAB(t), .2*t)
  )  
 
res = sym.dsolve(equations)

But I get this error:

  File "/usr/local/lib/python2.7/dist-packages/sympy/solvers/ode.py", line 1349, in classify_sysode
    t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0]
IndexError: list index out of range

Can anyone suggest where I am going wrong?  I think my syntax is not right.
Thanks.
Reply all
Reply to author
Forward
0 new messages