# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'

# import everything from SymPy.
from sympy import *

# Set up Jupyter notebook to display math.
init_printing() 

"""
Explanation: Modeling and Simulation in Python
Chapter 9
Copyright 2017 Allen Downey
License: Creative Commons Attribution 4.0 International
End of explanation
"""


from sympy.printing import latex

def show(expr, show_latex=False):
    """Display a SymPy expression.
    
    expr: SymPy expression
    show_latex: boolean
    """
    if show_latex:
        print(latex(expr))
    return expr

"""
Explanation: The following displays SymPy expressions and provides the option of showing results in LaTeX format.
End of explanation
"""


t = symbols('t')

"""
Explanation: Analysis with SymPy
Create a symbol for time.
End of explanation
"""


expr = t + 1

"""
Explanation: If you combine symbols and numbers, you get symbolic expressions.
End of explanation
"""


type(expr)

"""
Explanation: The result is an Add object, which just represents the sum without trying to compute it.
End of explanation
"""


expr.subs(t, 2)

"""
Explanation: subs can be used to replace a symbol with a number, which allows the addition to proceed.
End of explanation
"""


f = Function('f')

"""
Explanation: f is a special class of symbol that represents a function.
End of explanation
"""


type(f)

"""
Explanation: The type of f is UndefinedFunction
End of explanation
"""


f(t)

"""
Explanation: SymPy understands that f(t) means f evaluated at t, but it doesn't try to evaluate it yet.
End of explanation
"""


dfdt = diff(f(t), t)

type(dfdt)

"""
Explanation: diff returns a Derivative object that represents the time derivative of f
End of explanation
"""


alpha = symbols('alpha')

"""
Explanation: We need a symbol for alpha
End of explanation
"""


eq1 = Eq(dfdt, alpha*f(t))

"""
Explanation: Now we can write the differential equation for proportional growth.
End of explanation
"""


solution_eq = dsolve(eq1)

"""
Explanation: And use dsolve to solve it.  The result is the general solution.
End of explanation
"""


C1, p_0 = symbols('C1 p_0')

particular = solution_eq.subs(C1, p_0)

"""
Explanation: We can tell it's a general solution because it contains an unspecified constant, C1.
In this example, finding the particular solution is easy: we just replace C1 with p_0
End of explanation
"""


r, K = symbols('r K')

"""
Explanation: In the next example, we have to work a little harder to find the particular solution.
Solving the quadratic growth equation
We'll use the (r, K) parameterization, so we'll need two more symbols:
End of explanation
"""


eq2 = Eq(diff(f(t), t), r * f(t) * (1 - f(t)/K))

"""
Explanation: Now we can write the differential equation.
End of explanation
"""


solution_eq = dsolve(eq2)

"""
Explanation: And solve it.
End of explanation
"""


general = solution_eq.rhs

"""
Explanation: The result, solution_eq, contains rhs, which is the right-hand side of the solution.
End of explanation
"""


at_0 = general.subs(t, 0)

"""
Explanation: We can evaluate the right-hand side at $t=0$
End of explanation
"""


solutions = solve(Eq(at_0, p_0), C1)
type(solutions), len(solutions)

value_of_C1 = solutions[0]

"""
Explanation: Now we want to find the value of C1 that makes f(0) = p_0.
So we'll create the equation at_0 = p_0 and solve for C1.  Because this is just an algebraic identity, not a differential equation, we use solve, not dsolve.
The result from solve is a list of solutions.  In this case, we have reason to expect only one solution, but we still get a list, so we have to use the bracket operator, [0], to select the first one.
End of explanation
"""


particular = general.subs(C1, value_of_C1)

"""
Explanation: Now in the general solution, we want to replace C1 with the value of C1 we just figured out.
End of explanation
"""


particular = simplify(particular)

"""
Explanation: The result is complicated, but SymPy provides a method that tries to simplify it.
End of explanation
"""


particular.subs(t, 0)

"""
Explanation: Often simplicity is in the eye of the beholder, but that's about as simple as this expression gets.
Just to double-check, we can evaluate it at t=0 and confirm that we get p_0
End of explanation
"""


A = (K - p_0) / p_0

logistic = K / (1 + A * exp(-r*t))

"""
Explanation: This solution is called the logistic function.
In some places you'll see it written in a different form:
$f(t) = \frac{K}{1 + A e^{-rt}}$
where $A = (K - p_0) / p_0$.
We can use SymPy to confirm that these two forms are equivalent.  First we represent the alternative version of the logistic function:
End of explanation
"""


simplify(particular - logistic)

"""
Explanation: To see whether two expressions are equivalent, we can check whether their difference simplifies to 0.
End of explanation
"""


# Solution goes here

# Solution goes here

# Solution goes here

# Solution goes here

# Solution goes here

# Solution goes here

# Solution goes here

"""
Explanation: This test only works one way: if SymPy says the difference reduces to 0, the expressions are definitely equivalent (and not just numerically close).
But if SymPy can't find a way to simplify the result to 0, that doesn't necessarily mean there isn't one.  Testing whether two expressions are equivalent is a surprisingly hard problem; in fact, there is no algorithm that can solve it in general.
Exercises
Exercise: Solve the quadratic growth equation using the alternative parameterization
$\frac{df(t)}{dt} = \alpha f(t) + \beta f^2(t) $
End of explanation
"""


