I hope to take a few derivatives in time of a system of ordinary differential equations. So for instance with the following starting point:

from sympy import *

t = Symbol('t')
s = Function('s')(t)
i = Function('i')(t)
beta = Symbol('beta')
gamma = Symbol('gamma')

eqs = [Eq(s.diff(t), -beta*s*i), Eq(i.diff(t), beta*s*i-gamma*i)]

I'd like to get the second derivative with respect to time for s and for i, and higher time derivatives as well (as a function of s, i, beta, and gamma).

The natural symbolic solution seemed to me to be applying diff to the list or to the equations in the list. Neither of these appears to work, and after searching on stackoverflow and google I have no indications for what else I can even try.

Version info:

Sympy 1.10.1

Python 3.8.13

CodePudding user response：

As you have seen:

• it's not possible to apply a derivative to a Python list: diff doesn't know how to treat objects of type list, because they are not symbolic objects.
• it's not possible to apply a derivative to an Equality object (Eq), as these objects represents an equality relations and not a mathematical equation, so they do not support mathematical operations. :|

The best way to proceed is to write the equations in the form of (Left Hand Side) - (Right Hand Side):

eq1 = s.diff(t)   beta*s*i
eq2 = i.diff(t) - beta*s*i-gamma*i

Then, you can differentiate them. For example:

eq1.diff(t)
# output: beta*i(t)*Derivative(s(t), t)   beta*s(t)*Derivative(i(t), t)   Derivative(s(t), (t, 2))

CodePudding user response：

You already have the first derivatives as the rhs of your equations. You can just differentiate that again to get second derivatives:

In [10]: for eq in eqs:
    ...:     pprint(eq.rhs)
    ...: 
-β⋅i(t)⋅s(t)
β⋅i(t)⋅s(t) - γ⋅i(t)

In [11]: for eq in eqs:
    ...:     pprint(eq.rhs.diff(t))
    ...: 
         d                 d       
- β⋅i(t)⋅──(s(t)) - β⋅s(t)⋅──(i(t))
         dt                dt      
       d                 d            d       
β⋅i(t)⋅──(s(t))   β⋅s(t)⋅──(i(t)) - γ⋅──(i(t))
       dt                dt           dt  

If you want to use the differential equations to eliminate s'(t) and i'(t) from there then you can do that like:

In [12]: rep = {eq.lhs: eq.rhs for eq in eqs}

In [13]: for eq in eqs:
    ...:     pprint(eq.rhs.diff(t).subs(rep))
    ...: 
 2  2                                        
β ⋅i (t)⋅s(t) - β⋅(β⋅i(t)⋅s(t) - γ⋅i(t))⋅s(t)
   2  2                                                                   
- β ⋅i (t)⋅s(t)   β⋅(β⋅i(t)⋅s(t) - γ⋅i(t))⋅s(t) - γ⋅(β⋅i(t)⋅s(t) - γ⋅i(t))