Hi...
I'm very new to Sage, and I've run into some trouble. I have the following
system of equations, which come from a model of a simple protein
phosphorylation network:
s00,s01,s10,s11,k,p = var('s00 s01 s10 s11 k p')
eq1 = 0.55*k*s00 + 0.6*k*s01 + 0.6*k*s10 + 0.6*p*s01 + 0.6*p*s10 +
0.55*p*s11 + 33*s00 + 33*s01 + 33*s10 + 33*s11 == 33.0
eq2 = 0.55*k*s00 + 0.6*k*s01 + 0.6*k*s10 + 33*k == 0.33
eq3 = 0.6*p*s01 + 0.6*p*s10 + 0.55*p*s11 + 33*p == 0.33
eq4 = 0.11*k*s00 == 0.06*p*s01 + 0.06*p*s10
eq5 = (0.06*k + 0.06*p)*s01 == 0.055*k*s00 + 0.055*p*s11
eq6 = (0.06*k + 0.06*p)*s10 == 0.055*k*s00 + 0.055*p*s11
I solve it using:
slns =
solve([eq1,eq2,eq3,eq4,eq5,eq6],s00,s01,s10,s11,k,p,solution_dict=True)
and get back four solutions. My variables s00,s01,s10,s11,k,p are all
concentrations, and they often come back with negative or complex values; I
never find a solution where all are real and positive. I'd previously
integrated the underlying system of ODEs and found a real positive solution,
which I checked using:
print "\n\nFrom Runge-Kutta:"
intsol =
{s00:0.2608024056,s01:0.2390688718,s10:0.2390688718,s11:0.2608024056,k:0.0098712774,p:0.0098712774}
for x in [eq1,eq2,eq3,eq4,eq5,eq6]:
print x.subs(intsol).lhs().n(30),'==',x.subs(intsol).rhs().n(30)
and this is in fact a valid solution of the system of equations. To make
matters worse, I went back and checked the solutions that had been returned
by solve:
for s in slns:
print "\n\nCalculated solution #",slns.index(s)
print "s00:",s[s00].n(30),"\ts01:",s[s01].n(30),"\ts10:",s[s10].n(30)
print "s11:",s[s11].n(30),"\tk:",s[k].n(30),"\tp:",s[p].n(30)
for x in [eq1,eq2,eq3,eq4,eq5,eq6]:
print x.subs(s).lhs().n(30),'==',x.subs(s).rhs().n(30)
Only one of the solutions (solution #0) appears to be exactly correct; 2 and
3 are approximately correct, and 1 is, as far as I can tell, extremely
incorrect. So my questions are as follows:
1. Is solve() guaranteed to return all of the valid solutions for a system
like this one? What should I make of the fact that it seemed to miss one
(the one I care about, no less)?
2. Does the way that I'm double-checking solutions make sense? If so, then
why does one of the solutions look so wrong?
3. Am I using solve() correctly?
4. As an aside, I tried to restrict the domain of solutions by including:
for x in [s00,s01,s10,s11,k,p]:
assume(x>0)
assume(x,'real')
before the call to solve(). This appeared to make no difference at all with
regards to the solutions obtained, but it slowed my script down dramatically
-- from 3.2s without restricting the domain to 39.3s with the restrictions.
Is there a better way to specify that I only want real positive solutions,
or should I just throw out the ones I don't need after the call to solve()?
Thanks!
--craig
--
Craig Jolley
Foreign Postdoctoral Researcher
RIKEN Center for Developmental Biology
Laboratory for Systems Biology
http://www.cdb.riken.jp/lsb/index.html
--
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org
