CODEWITHSUNDEEP
pythonsympyode
SCLeo
SCLeo

Reputation: 383

Sympy not solving linear ODEs

I am trying to solve some very basic ODE systems:

t = symbols("t")
x, y = symbols("x y", cls=Function)
a = dsolve_system([
  Eq(y(t).diff(t), 1),
  Eq(x(t).diff(t), 1),
])

Output:

[[Eq(x(t), C1 + Integral(1, t)), Eq(y(t), C2 + Integral(1, t))]]

I am expecting something like x(t) = C1 + t and y(t) = C2 + t, but for some reason, the integrals are not getting solved.

When I tried to provide initial conditions:

t = symbols("t")
x, y = symbols("x y", cls=Function)
a = dsolve_system([
  Eq(y(t).diff(t), 1),
  Eq(x(t).diff(t), 1),
], ics={x(0): 1, y(0): 2})

The program just errored out:

NotAlgebraic                              Traceback (most recent call last)
<ipython-input-38-86edc5171650> in <module>()
      4   Eq(y(t).diff(t), 1),
      5   Eq(x(t).diff(t), 1),
----> 6 ], ics={x(0): 1, y(0): 2})
      7 a

14 frames
/usr/local/lib/python3.7/dist-packages/sympy/polys/numberfields.py in _minpoly_compose(ex, x, dom)
    583         res = _minpoly_rootof(ex, x)
    584     else:
--> 585         raise NotAlgebraic("%s doesn't seem to be an algebraic element" % ex)
    586     return res
    587 

NotAlgebraic: Integral(1, (t, 0)) doesn't seem to be an algebraic element

Upvotes: 0

Views: 285

Answers (2)

Oscar Benjamin
Oscar Benjamin

Reputation: 14480

One difference between dsolve_system and dsolve (which internally calls dsolve_system) is that in dsolve_system integrals are not automatically computed unless you pass doit=True:

In [16]: dsolve_system(eqs)
Out[16]: 
⎡⎡            ⌠                   ⌠     ⎤⎤
⎢⎢y(t) = C₁ + ⎮ 1 dt, x(t) = C₂ + ⎮ 1 dt⎥⎥
⎣⎣            ⌡                   ⌡     ⎦⎦

In [17]: dsolve_system(eqs, doit=True)
Out[17]: [[y(t) = C₁ + t, x(t) = C₂ + t]]

In this case it doesn't make much difference but in general computing the integrals can be a lot slower than the rest of the calculations so it is not done by default when calling dsolve_system.

You can also call doit yourself to evaluate the integrals:

In [18]: [[sol1, sol2]] = dsolve_system(eqs)

In [19]: sol1
Out[19]: 
            ⌠     
y(t) = C₁ + ⎮ 1 dt
            ⌡     

In [20]: sol1.doit()
Out[20]: y(t) = C₁ + t

The problem with dsolve_system and the ics argument is a bug to do with how the unevaluated integrals get handled internally. The bug is not seen with dsolve because it always evaluates the integrals.

Upvotes: 1

SCLeo
SCLeo

Reputation: 383

For anyone who encountered this issue in the future, here is how I workaround it:

Instead of using dsolve_system, I used dsolve. For whatever reason, dsolve does not produce those weird Integral(1,t)'s. However, providing initial conditions will still error out. Therefore, I manually solved C1 and C2 using solve. I am not sure how robust my methods are, but at least they are working:

t = symbols("t")
x, y = symbols("x y", cls=Function)

solutions = dsolve([
  Eq(x(t).diff(t), 1),
  Eq(y(t).diff(t), 1),
])
C1, C2 = symbols("C1 C2")
initial_conditions = solve([
  Eq(3, solutions[0].rhs.subs(t, 0)), # 3 == the first function evaluated at t = 0
  Eq(5, solutions[1].rhs.subs(t, 0)), # 5 == the second function evaluated at t = 0
])
solutions = [
  solutions[0].subs(C1, initial_conditions[C1]).subs(C2, initial_conditions[C2]),
  solutions[1].subs(C1, initial_conditions[C1]).subs(C2, initial_conditions[C2]),
]
solutions

Output:

[Eq(x(t), t + 3), Eq(y(t), t + 5)]

Nevertheless, if you know what I did wrong originally or why doesn't providing initial conditions work when calling dsolve, please let me know.

Upvotes: 1

Related Questions