How to solve this logistic differential in sage?
y=function('y')(x)
(b,r,y0)=(30,0.8,1.5)
de=diff(y,x)==r*y*((1-y)/b) # logistic equation
soln=desolve(de,y,ics=[0,y0])
mathematica solution is: y(x)=(b*exp(r*x)*y0)/(exp(r*x)*y0-y0+b)
Other questions/answers found on the logistic topic were not helpful for this problem
The posted code:
( b, r, y0 ) = ( 30, 8/10, 3/2 )
var( 'x' );
y = function('y')(x)
de = diff( y, x ) == r*y*(1-y)/b
sol = desolve( de, y, ics=[0,y0] )
already leads to the solution (using implicit functions of y and x):
sage: sol
-75/2*log(y(x) - 1) + 75/2*log(y(x)) == x + 75/2*log(2) + 75/2*log(3/2)
In the ode-world this is already a solution. (Although i had a similar situation in an exam, and the examiner considered
i stopped too early, since the equation can be solved explicitly.) Here, the coefficients of the two log expressions involving a non-constant function of y(x) are opposite, $\pm 75/2$, so using exp on both sides the expression $(y-1)/y=1-1/y$ leads to an algebraic equation of first order in $y$, that can be solved. If we want to do this in sage, there is one more step needed, before calling solve (for a combined exponential/logarithmic and algebraic equation), we substitute with an other variable, z, the expression y(x). The adjusted code for this plan and a possible answer to the question is as follows:
( b, r, y0 ) = ( 30, 8/10, 3/2 )
var( 'x,z' );
y = function('y')(x)
de = diff( y, x ) == r*y*(1-y)/b
sol = desolve( de, y, ics=[0,y0] )
sol = sol . subs( { y(x): z } )
Y = sol . solve( z, to_poly_solve=True )[0].rhs()
print "Solution: Y=%s" % Y
print "Check is", bool( diff(Y,x) == r*Y*(1-Y)/ b )
print "Y(0) = %s" % Y.subs( {x:0} )
Results:
Solution: Y=3*e^(2/75*x)/(3*e^(2/75*x) - 1)
Check is True
Y(0) = 3/2
Related question, answers and discussion: logistic differential equation
The above formula corresponds to...
sage: ( b, r, y0 ) = ( 30, 8/10, 3/2 )
sage: Y = b*exp(r*x)*y0/(exp(r*x)*y0-y0+b)
sage: Y
30*e^(4/5*x)/(e^(4/5*x) + 19)
But please observe the difference between the two expressions:
de=diff(y,x)==r*y*((1-y)/b) as posted, and
de=diff(y,x)==r*y*(1-y/b) as leading to the above solution. Let's see if the code fits in place:
( b, r, y0 ) = ( 30, 8/10, 3/2 )
var( 'x,z' );
y = function('y')(x)
de = diff( y, x ) == r*y*(1-y/b)
sol = desolve( de, y, ics=[0,y0] )
Y = sol . subs( { y(x): z } ) . solve( z, to_poly_solve=True )[0].rhs()
print "Solution: Y=%s" % Y
which gives indeed:
Solution: Y=30*e^(4/5*x)/(e^(4/5*x) + 19)
dan_fulea: Thank you. Solved my problem/question. I am a senior (very) not a student. Math is my hobby. My mathematica version does not run on a 64 bit machine and therefore, switching to sage and learning on my own with the help of a few texts books. I have not seen/learned the code "sol.subs({y(x):z}).solve(z,to_poly_solve=True)[0].rhs() in any of my text books. After your help I found with 'solve? evaluate' but no explanation as to why. A bit exotic for my level. The code fixes my problem. Thanks again. כול הכבוד
P.S. I miss-typed my original questions. Meant to write (1-y/b); not (1-y)/b. That error caused some unecessary confusion.
Asked: 2017-09-15 21:37:47 +0100
Seen: 595 times
Last updated: Sep 16 '17
Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.
