problem when using dsolve() to solve ordinary differential equations

[itstyren]

When I use dsolve() like this, I get an error in the last
How can I solve it, thank you！

import sympy as sy

u1_t = 0  # u_0=0
u2_t = 0  # u_0=0

def differential_equation_x1(t, x_1, x_2, u1_t, u2_t):
    return sy.diff(x_1(t), t,
                   1) + x_1(t) + x_2(t) - u1_t - u2_t  


def differential_equation_x2(t, x_1, x_2, u1_t, u2_t):
    return sy.diff(x_2(t), t,
                   1) - x_1(t) - x_2(t) + u1_t - u2_t 


t = sy.symbols('t')  
x_1 = sy.Function('x_1')  
x_2 = sy.Function('x_2')  
x_t = sy.dsolve((differential_equation_x1(t, x_1, x_2, u1_t, u2_t),
                 differential_equation_x2(t, x_1, x_2, u1_t, u2_t)),
                (x_1(t), x_2(t)),
                ics={
                    x_1(0): 10,
                    x_2(0): 5
                })
sy.pprint(x_t)


def differential_equation_l1(t, l_1, l_2, x_t):
    return sy.diff(l_1(t), t, 1) + x_t[0].rhs - l_1(t)+l_2(t)

def differential_equation_l2(t, l_1, l_2, x_t):
    return sy.diff(l_2(t), t, 1) + x_t[1].rhs - l_1(t) + l_2(t)

t = sy.symbols('t')  
l_1 = sy.Function('l_1')  
l_2 = sy.Function('l_2')  

l_t = sy.dsolve((differential_equation_l1(t, l_1, l_2, x_t),
                 differential_equation_l2(t, l_1, l_2, x_t)), (l_1(t), l_2(t)),
                ics={
                    l_1(1): 0,
                    l_2(1): 0
                },
                )
sy.pprint(l_t)

[x₁(t) = -15⋅t + 10, x₂(t) = 15⋅t + 5]
Traceback (most recent call last):
File "example_3.py", line 77, in
l_2(1): 0
File "D:\Anaconda3\lib\site-packages\sympy\solvers\ode.py", line 609, in dsolve
raise NotImplementedError
NotImplementedError

[oscarbenjamin]

Thanks for reporting this. Here is simpler code to reproduce this:

from sympy import *

t = symbols('t')
l_1 = Function('l_1')
l_2 = Function('l_2')
eqs = [
        Eq(l_1(t).diff(t), + l_1(t) - l_2(t) + 15*t - 10),
        Eq(l_2(t).diff(t), + l_1(t) - l_2(t) - 15*t - 5),
    ]

pprint(eqs[0])
pprint(eqs[1])

l_t = dsolve(eqs, [l_1(t), l_2(t)])
pprint(l_t)

The ODEs are:

d                                    
──(l₁(t)) = 15⋅t + l₁(t) - l₂(t) - 10
dt                                   
d                                    
──(l₂(t)) = -15⋅t + l₁(t) - l₂(t) - 5
dt  

Traceback:

Traceback (most recent call last):
  File "myfiles/g.py", line 14, in <module>
    l_t = dsolve(eqs, [l_1(t), l_2(t)])
  File "/Users/enojb/current/sympy/sympy/sympy/solvers/ode.py", line 619, in dsolve
    raise NotImplementedError
NotImplementedError

These should be solvable as the system is constant coefficient non-homogeneous (CCNH). Various special cases of CCNH are solvable for systems of 2 or 3 equations. These can be solved in general using the matrix exponential but it isn't implemented yet in dsolve.

See here:
https://en.wikipedia.org/wiki/Matrix_exponential#Inhomogeneous_case_generalization:_variation_of_parameters

[bgoodman44]

I'm having the same issue on the following system of ODEs:

import sympy
sympy.init_printing()

t = sympy.Symbol('t',real=True)
x = sympy.Function('x')(t)
v = sympy.Function('v')(t)
w = sympy.Symbol('omega',real=True,positive=True,nonzero=True)

eq1 = sympy.Eq(x.diff(),v)
eq2 = sympy.Eq(v.diff(),sympy.cos(w*t))
eqs = [eq1,eq2]

sympy.dsolve(eqs,[x,v])

System pprint:

d              
──(x(t)) = v(t)
dt 

d                  
──(v(t)) = cos(ω⋅t)
dt    

[oscarbenjamin]

That's a different problem actually. The dsolve system code lacks the capability to solve an integrable system of ODEs. It needs to the equivalent of the nth_algebraic method but for systems.

[bgoodman44]

Ouch. That basically eliminates most/many/all physics based ODEs describing dynamic systems.

Was hoping to use sympy as a replacement for Matlab symbolic math, but doesn't sound like the maturity is there just yet.

[oscarbenjamin]

You've misunderstood what I meant by "integrable" I think. What I mean is: trivial ODE systems that can be solved simply by integration. There isn't currently a solver implemented just for that. It wouldn't be hard to add though.

[bgoodman44]

Nah, we're on the same page, and I understand the work-around. I'd imagine this would be a common use-case though, and would be a useful feature.

[oscarbenjamin]

If anyone wants to write a pull request this we just need to do the same as nth_algebraic does for single ODEs. The nth_algebraic solver replaces derivatives with a "function" whose inverse is an integral:

sympy/sympy/solvers/ode.py

Lines 4574 to 4579 in ac1eaa2

	class diffx(Function):
	def inverse(self):
	# don't use integrate here because fx has been replaced by _t
	# in the equation; integrals will not be correct while solve
	# is at work.
	return lambda expr: Integral(expr, var) + Dummy('C')

Then the system of ODEs can be solved by solve as if it were a system of algebraic equations:

sympy/sympy/solvers/ode.py

Line 4623 in ac1eaa2

solns = solve(subs_eqn, func, simplify=False)

We already have that implemented for single ODEs so the code there can be reused. We just need to find a place in the systems matching code to try this and see if it works.

[mijo2]

@oscarbenjamin I missed this issue when we were working in the type 2 solver. This issue can be closed now.

In [232]: t = symbols('t') 
     ...: l_1 = Function('l_1') 
     ...: l_2 = Function('l_2') 
     ...: eqs = [ 
     ...:         Eq(l_1(t).diff(t), + l_1(t) - l_2(t) + 15*t - 10), 
     ...:         Eq(l_2(t).diff(t), + l_1(t) - l_2(t) - 15*t - 5), 
     ...:     ]                                                                                                                                                                                             

In [233]: dsolve(eqs)                                                                                                                                                                                       
Out[233]: 
⎡                                                               ⌠                                                                  ⌠                        ⎤
⎢                           ⌠                 ⌠                 ⎮ ⎛      2           ⎞                           ⌠                 ⎮ ⎛      2           ⎞   ⎥
⎢l₁(t) = C₁ + C₂⋅t + C₂ + t⋅⎮ (30⋅t - 5) dt + ⎮ (30⋅t - 5) dt + ⎮ ⎝- 30⋅t  - 10⋅t - 5⎠ dt, l₂(t) = C₁ + C₂⋅t + t⋅⎮ (30⋅t - 5) dt + ⎮ ⎝- 30⋅t  - 10⋅t - 5⎠ dt⎥
⎣                           ⌡                 ⌡                 ⌡                                                ⌡                 ⌡                        ⎦

In [234]: neq_nth_linear_constant_coeff_match(eqs, [l_1(t), l_2(t)], t)                                                                                                                                     
Out[234]: 
{'no_of_equation': 2,
 'eq': [Eq(Derivative(l_1(t), t), 15*t + l_1(t) - l_2(t) - 10),
  Eq(Derivative(l_2(t), t), -15*t + l_1(t) - l_2(t) - 5)],
 'func': [l_1(t), l_2(t)],
 'order': {l_1(t): 1, l_2(t): 1},
 'is_linear': True,
 'is_constant': True,
 'is_homogeneous': False,
 'is_general': True,
 'func_coeff': Matrix([
 [-1, 1],
 [-1, 1]]),
 'rhs': Matrix([
 [15*t - 10],
 [-15*t - 5]]),
 'type_of_equation': 'type2'}

In [235]: sol = dsolve(eqs)                                                                                                                                                                                 

In [236]: checksysodesol(eqs, sol)                                                                                                                                                                          
Out[236]: (True, [0, 0])

[mijo2]

I will add a regression test case for the same in the current PR.

[oscarbenjamin]

It can be closed when the regression test is merged to master