# solution to logistic differential y'=r*y*(1-y)/b

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

edit retag close merge delete

Sort by » oldest newest most voted

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 ).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

more

re:logistic function If I understand your reply correctly, y(x)=3e^(2/75x)/(3e^(2/75x) - 1) But y(x)=3e^(2/75x)/(3e^(2/75x) - 1) is not a logistic function Plot it out and plot out y(x)=(bexp(rx)y0)/(exp(rx)*y0-y0+b) Again thank you for responding.

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 ).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).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.