# Sympy Computation with dsolve for Logistic Differential Equation

**URL:** <https://discourse.julialang.org/t/sympy-computation-with-dsolve-for-logistic-differential-equation/104338>\
**Category:** General Usage\
**Tags:** question\
**Created:** [September 28, 2023, 7:48am UTC](https://discourse.julialang.org/t/sympy-computation-with-dsolve-for-logistic-differential-equation/104338 "2023-09-28T07:48:02Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Freya\_the\_Goddess](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/freya_the_goddess/32/36835_2.png) [@Freya\_the\_Goddess](https://discourse.julialang.org/u/Freya_the_Goddess)\
**Post date:** [September 28, 2023, 7:48am UTC](https://discourse.julialang.org/t/sympy-computation-with-dsolve-for-logistic-differential-equation/104338/1 "2023-09-28T07:48:02Z")

</div>

I was computing the solution for:  
y' = ky(L-y)

the solution should be:

y(t) = C\_{1} e^{0.6t} (2000 - y(t))

the value of C\_{1} should be 2/3. Thus

y(t) = \frac{4000/3}{2/3 + e^{-0.6t}}

I believe I have followed the formula correctly and enter into the code, but it seems to be all wrong, this is the code:

```julia

using SymPy, Plots, LaTeXStrings, Plots.PlotMeasures
@syms q, t, r
@syms q()

# Computing IVP with SymPy
q(t).diff(t) 
k = 0.0003
L = 2000
# Let Q' = kQ (L - Q)
diffeq = Eq(q(t).diff(t), k*q(t)*(L-q(t)) ); string(diffeq)
# To solve the ODE, pass it and the function to solve for to dsolve.
qt = dsolve(diffeq, q(t))
# To solve the Initial Value Problem q(0)=800
ivp = dsolve(diffeq, q(t), ics = Dict(q(0) => 800))

plot(y, 0,10,
	legend=:topright, bottom_margin=5mm, 
	left_margin=5mm,
	xlabel="years", ylabel="population",
	label=L"y(t) = \frac{4000/3}{2/3 + e^{0.6 t}}",
	size=(720, 360), tickfontsize=10)

```

![Screenshot from 2023-09-28 07-39-07](https://global.discourse-cdn.com/julialang/original/3X/5/8/5826d33f8ab89a5792a37d0b32f05ed4a6798076.png)

What did I do wrong or is it because Logistic Differential Equation can’t be solved with dsolve?

---

<div class="post-metadata">

**Author:** ![j\_verzani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_verzani/32/8551_2.png) [@j\_verzani](https://discourse.julialang.org/u/j_verzani)\
**Post date:** [September 28, 2023, 11:20am UTC](https://discourse.julialang.org/t/sympy-computation-with-dsolve-for-logistic-differential-equation/104338/2 "2023-09-28T11:20:28Z")

</div>

You might try it this way:

```julia
julia> @syms K, L, y(), t
(K, L, y, t)

julia> ∂ = Differential(t)
Differential(t)

julia> u = dsolve(∂(y(t)) ~ K*y(t)*(L - y(t)))
           L⋅(C₁ + K⋅t) 
        L⋅ℯ             
y(t) = ─────────────────
        L⋅(C₁ + K⋅t)    
       ℯ - 1

julia> ex =subs(rhs(u), L=>2000, K => 3//10000)
                3⋅t
      2000⋅C₁ + ───
                 5 
2000⋅ℯ             
───────────────────
            3⋅t    
  2000⋅C₁ + ───    
             5     
 ℯ - 1

```

If you have an initial condition for `y(0)` you can pass it to `dsolve`. If you want to substitute in for `C_1` directly, you need to get that symbol. The command `first(free_symbols(ex))` identifies it.

---

<div class="post-metadata">

**Author:** ![Freya\_the\_Goddess](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/freya_the_goddess/32/36835_2.png) [@Freya\_the\_Goddess](https://discourse.julialang.org/u/Freya_the_Goddess)\
**Post date:** [September 30, 2023, 1:09pm UTC](https://discourse.julialang.org/t/sympy-computation-with-dsolve-for-logistic-differential-equation/104338/3 "2023-09-30T13:09:16Z")

</div>

I tried till this:

```julia

using SymPy, Plots, LaTeXStrings, Plots.PlotMeasures
@syms K, L, y(), t

∂ = Differential(t)
yt = dsolve(∂(y(t)) ~ K*y(t)*(L - y(t)))
sol_yt = subs(rhs(yt), L=>2000, K => 3//10000)

# Manual way
# Find the constant C₁
c1 = free_symbols(yt)[1]
# To solve the Initial Value Problem, initial condition y(0) = 800
C1 = solve(sol_yt(t=>0) ~ 800, c1)
#ivp_yt = sol_yt(c1=>C1[1])

```

the computation to get constant C\_{1} is taking too long.

But I made it with Python code, even if it is not Julia related here

```julia
# https://natelemoine.com/using-sympy-for-biological-odes/

import sympy as sm
import numpy as np
import matplotlib.pyplot as plt

from sympy.abc import L, k, t
y = sm.Function('y')(t)

# The differential equation: y' = ky(L-y)
# Define the derivative of y with respect to t (the left-hand side of the ODE),
# and then define the right-hand side of the ODE.
dy = y.diff(t)
rhs = k*y*(L - y)

# Define the equality
eq = sm.Eq(dy, rhs)
print('the Differential equation is:')
sm.pretty_print(eq)

print('the Differential equation solution is:')
sol = sm.dsolve(eq)
sm.pretty_print(sol)

print('the solution of the Differential equation for initial value of y(0) = 800:')

t0 = sol.args[1].subs({'t': 0})
n0 = 800
eq_init = sm.Eq(n0, t0)

C1 = t0.args[2].args[0].args[0]
t0_sol = sm.solve(eq_init, C1)
# substitute the expression for C1 back into the equation
final = sol.args[1].subs(C1, t0_sol[0])
final_eq = final.simplify()
sm.pretty_print(final_eq.subs([(L,2000),(k,0.0003)]))

f_num = sm.lambdify([k, L, t], final_eq)
t_coords = np.linspace(0, 10, 100)
y_coords = f_num(0.0003, 2000, t_coords)

# Plot with SymPy Plotting Backends module:
sm.plot(final_eq.subs([(L,2000),(k,0.0003)]), (t, 0, 10),
	xlabel='$t$ (year)', ylabel='$y(t)$ (population)', title='Logistic Differential Equation Plot')

# Plot with matplotlib
plt.figure()
plt.plot(t_coords, y_coords)
plt.title('Logistic Differential Equation Plot')
plt.xlabel("$t$ (year)")
plt.ylabel("$y(t)$ (population)")

plt.show()

```

---

<div class="post-metadata">

**Author:** ![j\_verzani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_verzani/32/8551_2.png) [@j\_verzani](https://discourse.julialang.org/u/j_verzani)\
**Post date:** [September 30, 2023, 4:02pm UTC](https://discourse.julialang.org/t/sympy-computation-with-dsolve-for-logistic-differential-equation/104338/4 "2023-09-30T16:02:32Z")

</div>

If you know the initial value, you should let SymPy do the work within dsolve:

```julia
using SymPy
@syms K, L, y(), t
∂ = Differential(t)
yt = dsolve(∂(y(t)) ~ K*y(t)*(L - y(t)), ics=Dict(y(0) => 800)) # equation of t; K, L
yt = rhs(yt) # work with just the expression, not the equation
yt = subs(yt, L=>2000, K => 3//10000) # leaves as an expression of t alone

```

The value of `yt` can be passed on to `plot` using the recipe, or after calling `lambdify` can be passed as a function to other actions.

The `ics` argument allows you to pass a dictionary, as shown.

The earlier answer suggested what to do if you knew the value of `C1` (find the variable, the substitute in for the known value). That isn’t suggested unless that is all you have to work with.

---

<div class="post-metadata">

**Author:** ![Freya\_the\_Goddess](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/freya_the_goddess/32/36835_2.png) [@Freya\_the\_Goddess](https://discourse.julialang.org/u/Freya_the_Goddess)\
**Post date:** [October 2, 2023, 11:53am UTC](https://discourse.julialang.org/t/sympy-computation-with-dsolve-for-logistic-differential-equation/104338/5 "2023-10-02T11:53:31Z")

</div>

I’m finally able to make it work with Julia code, thanks a lot

```julia
# https://discourse.julialang.org/t/sympy-computation-with-dsolve-for-logistic-differential-equation/104338/3

using SymPy, Plots, LaTeXStrings, Plots.PlotMeasures
@syms K, L, y(), t, y0

∂ = Differential(t)
yt_sym = dsolve(∂(y(t)) ~ K*y(t)*(L - y(t)), ics=Dict(y(0) => y0)) # equation of t; K, y0, L
yt_ivp = dsolve(∂(y(t)) ~ K*y(t)*(L - y(t)), ics=Dict(y(0) => 800)) # equation of t; K, L

yt = rhs(yt_ivp) # work with just the expression, not the equation
yt_final = subs(yt, L=>2000, K => 3//10000) # leaves as an expression of t alone

println("Logistic Differential Equation with Julia")
println("The solution")
sympy.pretty_print(yt_sym)

println("The solution for initial value problem y(0)=800")
sympy.pretty_print(yt_ivp)

println("Logistic Differential Equation with L=2000 and K = 0.00003")
sympy.pretty_print(yt_final)

plot(yt_final, 0,10,
	legend=:bottomright, bottom_margin=5mm, 
	left_margin=5mm,
	xlabel="years", ylabel="population",
	label=L"y(t) = \frac{4000/3}{2/3 + e^{0.6 t}}",
	size=(720, 360), tickfontsize=10)

```

 ![Capture d’écran_2023-10-02_18-50-40](https://global.discourse-cdn.com/julialang/original/3X/c/e/ce2f7814877887bbd593c39403ca12c44c0e3589.png)