KeyError in dsolve matching system of 1st order nonlinear ODEs

[oscarbenjamin]

This example comes from stackoverflow:
https://stackoverflow.com/questions/16909779/any-way-to-solve-a-system-of-coupled-differential-equations-in-python

The system of ODEs is

v11, v22, v12 = symbols('v11, v22, v12', cls=Function)
t = Symbol('t')
eq1 = Eq(v11(t).diff(t), -12*v12(t)**2)
eq2 = Eq(v22(t).diff(t), +12*v12(t)**2)
eq3 = Eq(v12(t).diff(t), 6*v11(t)*v12(t) - 6*v12(t)*v22(t) -36*v12(t))
dsolve([eq1, eq2, eq3], [v11(t), v22(t), v12(t)])

and the result:

Traceback (most recent call last):
  File "j.py", line 7, in <module>
    dsolve([eq1, eq2, eq3], [v11(t), v22(t), v12(t)])
  File "/Users/enojb/current/sympy/sympy/sympy/solvers/ode.py", line 595, in dsolve
    match = classify_sysode(eq, func)
  File "/Users/enojb/current/sympy/sympy/sympy/solvers/ode.py", line 1657, in classify_sysode
    type_of_equation = check_nonlinear_3eq_order1(eq, funcs, func_coef)
  File "/Users/enojb/current/sympy/sympy/sympy/solvers/ode.py", line 2015, in check_nonlinear_3eq_order1
    if eq[1].has(r1[F2]) and not eq[1].has(r1[F3]):
KeyError: F3_

Working on this a bit I get another error:

In [29]: repr([eq1, eq33])                                                                                                        
Out[29]: '[Eq(Derivative(v11(t), t), -12*v12(t)**2), Eq(Derivative(v12(t), t), -6*(C - 2*v11(t) + 6)*v12(t))]'

In [30]: dsolve([eq1, eq33], [v11(t), v12(t)])                                                                                    
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-30-6c88cb48f966> in <module>
----> 1 dsolve([eq1, eq33], [v11(t), v12(t)])

~/current/sympy/sympy/sympy/solvers/ode.py in dsolve(eq, func, hint, simplify, ics, xi, eta, x0, n, **kwargs)
    593     """
    594     if iterable(eq):
--> 595         match = classify_sysode(eq, func)
    596         eq = match['eq']
    597         order = match['order']

~/current/sympy/sympy/sympy/solvers/ode.py in classify_sysode(eq, funcs, **kwargs)
   1650             if matching_hints['no_of_equation'] == 2:
   1651                 if order_eq == 1:
-> 1652                     type_of_equation = check_nonlinear_2eq_order1(eq, funcs, func_coef)
   1653                 else:
   1654                     type_of_equation = None

~/current/sympy/sympy/sympy/solvers/ode.py in check_nonlinear_2eq_order1(eq, func, func_coef)
   1957     num, den = (
   1958         (r1[f].subs(x(t),u).subs(y(t),v))/
-> 1959         (r2[g].subs(x(t),u).subs(y(t),v))).as_numer_denom()
   1960     R1 = num.match(f1*g1)
   1961     R2 = den.match(f2*g2)

TypeError: 'NoneType' object is not subscriptable

[Awkwafina]

Hi! I am newbie in GSOC. So, I would like to work on this issue. I am a math major student and have some experience with Python. I hope I can handle it. Let me know if it still open.

[Yhozen]

I'm getting the same error with a basic SIR model. Any ideas?

[oscarbenjamin]

Can you show the equations that you are using?

Note that "basic" equations like this very often don't have any analytic solution so although there is a bug here there isn't necessarily any reason to expect a useful output.

[Yhozen]

That might be the reason, this is my code:

from sympy import *
from sympy.abc import t

# Contact rate, beta, and mean recovery rate, gamma, (in 1/days).
beta, gamma = 0.2, 1./10 

S = Function('S')(t)
I = Function('I')(t)
R = Function('R')(t)
 
Sdt = S.diff(t)
Idt = I.diff(t)
Rdt = R.diff(t)
 
Sdt_right =  -1 * beta * S * I 
Idt_right =  beta * S * I  - gamma * I 
Rdt_right = gamma * I
 

 
eq1 = Eq(Sdt, Sdt_right)
eq2 = Eq(Idt, Idt_right)
eq3 = Eq(Rdt, Rdt_right)
 
ODEs = [eq1,eq2,eq3]
 
dsolve(ODEs, [S,I,R])

[oscarbenjamin]

The ODEs are

In [21]: eq1                                                                                                                      
Out[21]: 
d                        
──(S(t)) = -0.2⋅I(t)⋅S(t)
dt                       

In [22]: eq2                                                                                                                      
Out[22]: 
d                                  
──(I(t)) = 0.2⋅I(t)⋅S(t) - 0.1⋅I(t)
dt                                 

In [23]: eq3                                                                                                                      
Out[23]: 
d                  
──(R(t)) = 0.1⋅I(t)
dt  

There are parametric solutions given here but they seem fairly complex:
https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#Exact_analytical_solutions_to_the_SIR_model
There is no algorithm in sympy that can solve this kind of system currently. The way to solve it might be:

1. Ignore eq3 since if you can solve the other equations for I then R will just be the integral of that.

2. Divide equations 1 and 2 and eliminate t from the autonomous system:

In [35]: s = Function('s')                                                                                                        

In [36]: i = Symbol('i')                                                                                                          

In [37]: eq4 = Eq(s(i).diff(i), (eq1.rhs/eq2.rhs).subs(S, s(i)).subs(I, i))                                                       

In [38]: eq4                                                                                                                      
Out[38]: 
d             -0.2⋅i⋅s(i)    
──(s(i)) = ──────────────────
di         0.2⋅i⋅s(i) - 0.1⋅i

3. Solve that for S(I):

In [39]: dsolve(eq4)                                                                                                              
Out[39]: 
             ⎛         2.0⋅i⎞
s(i) = -0.5⋅W⎝-2.0⋅C₁⋅ℯ     ⎠

4. Use this to eliminate either S or I from one of the original ODEs and then solve:

In [46]: sol = dsolve(nsimplify(eq4))                                                                                             

In [47]: sol                                                                                                                      
Out[47]: 
         ⎛    2⋅i⎞ 
       -W⎝C₁⋅ℯ   ⎠ 
s(i) = ────────────
            2      

In [48]: solve(sol.subs(s(i), S), i)                                                                                              
Out[48]: 
⎡                                          ⎛         -2⋅S(t) ⎞⎤
⎢   ⎛      _________________⎞              ⎜-2⋅S(t)⋅ℯ        ⎟⎥
⎢   ⎜     ╱        -2⋅S(t)  ⎟           log⎜─────────────────⎟⎥
⎢   ⎜    ╱  -S(t)⋅ℯ         ⎟   log(2)     ⎝        C₁       ⎠⎥
⎢log⎜-  ╱   ─────────────── ⎟ + ──────, ──────────────────────⎥
⎣   ⎝ ╲╱           C₁       ⎠     2               2           ⎦

In [49]: eq1p = eq1.subs(I, _[1])                                                                                                 

In [50]: eq1p                                                                                                                     
Out[50]: 
                        ⎛         -2⋅S(t) ⎞
d                       ⎜-2⋅S(t)⋅ℯ        ⎟
──(S(t)) = -0.1⋅S(t)⋅log⎜─────────────────⎟
dt                      ⎝        C₁       ⎠

In [51]: dsolve(eq1p)                                                                                                             
Out[51]: 
S(t)                                            
 ⌠                                              
 ⎮                1                             
 ⎮   ─────────────────────────── dy = C₂ - 0.1⋅t
 ⎮     ⎛   ⎛    -2⋅y ⎞         ⎞                
 ⎮     ⎜   ⎜-y⋅ℯ     ⎟         ⎟                
 ⎮   y⋅⎜log⎜─────────⎟ + log(2)⎟                
 ⎮     ⎝   ⎝    C₁   ⎠         ⎠                
 ⌡ 

That just ends up with the result given implicitly which is presumably similar to the parametric solutions from wikipedia.