New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

BUG: dsolve fails for linear first order ODE with three equations #15407

Open
Tillsten opened this Issue Oct 19, 2018 · 3 comments

Comments

Projects
None yet
4 participants
@Tillsten

Tillsten commented Oct 19, 2018

The following code returns False for c(t).

import sympy
from sympy import Eq, Derivative
a, b, c = sympy.symbols('a b c', cls=sympy.Function)
t, a_b, a_c = sympy.symbols('t a_b a_c', real=True)

eqs = [Eq(Derivative(a(t), t), (-a_b - a_c)*a(t)), Eq(Derivative(b(t), t), a_b*a(t)), Eq(Derivative(c(t), t), a_c*a(t))]
sympy.dsolve(eqs)

@Tillsten Tillsten changed the title BUG: Solve fails for linear first order ODE with three equations BUG: dsolve fails for linear first order ODE with three equations Oct 19, 2018

@Tillsten Tillsten closed this Oct 19, 2018

@Tillsten Tillsten reopened this Oct 19, 2018

skirpichev added a commit to skirpichev/diofant that referenced this issue Oct 19, 2018

@asmeurer

This comment has been minimized.

Member

asmeurer commented Oct 19, 2018

It looks like the code is not properly checking for a case where it divides by a quantity and it divides by zero, producing nan, and Eq(c(t), nan) reduces to False.

The issue is in _linear_3eq_order1_type1.

skirpichev added a commit to skirpichev/diofant that referenced this issue Oct 20, 2018

skirpichev added a commit to skirpichev/diofant that referenced this issue Oct 21, 2018

@oscarbenjamin

This comment has been minimized.

Contributor

oscarbenjamin commented Nov 5, 2018

The equations solved by _linear_3eq_order1_type1 are also solved by _linear_neq_order1_type1. This issue can be fixed by removing _linear_3eq_order1_type1.

I tried this:

diff --git a/sympy/solvers/ode.py b/sympy/solvers/ode.py
index 8fc6f51..5d54ad8 100644
--- a/sympy/solvers/ode.py
+++ b/sympy/solvers/ode.py
@@ -7828,7 +7828,7 @@ def sysode_linear_3eq_order1(match_):
             if not j.has(x(t), y(t), z(t)):
                 raise NotImplementedError("Only homogeneous problems are supported, non-homogenous are n
     if match_['type_of_equation'] == 'type1':
-        sol = _linear_3eq_order1_type1(x, y, z, t, r, eq)
+        sol = _linear_neq_order1_type1(match_)
     if match_['type_of_equation'] == 'type2':
         sol = _linear_3eq_order1_type2(x, y, z, t, r, eq)
     if match_['type_of_equation'] == 'type3':
@@ -7839,35 +7839,6 @@ def sysode_linear_3eq_order1(match_):
         sol = _linear_neq_order1_type1(match_)
     return sol
 
-def _linear_3eq_order1_type1(x, y, z, t, r, eq):
-    r"""
-    .. math:: x' = ax
-
-    .. math:: y' = bx + cy
-
-    .. math:: z' = dx + ky + pz
-
-    Solution of such equations are forward substitution. Solving first equations
-    gives the value of `x`, substituting it in second and third equation and
-    solving second equation gives `y` and similarly substituting `y` in third
-    equation give `z`.
-
-    .. math:: x = C_1 e^{at}
-
-    .. math:: y = \frac{b C_1}{a - c} e^{at} + C_2 e^{ct}
-
-    .. math:: z = \frac{C_1}{a - p} (d + \frac{bk}{a - c}) e^{at} + \frac{k C_2}{c - p} e^{ct} + C_3 e^{
-
-    where `C_1, C_2` and `C_3` are arbitrary constants.
-
-    """
-    C1, C2, C3, C4 = get_numbered_constants(eq, num=4)
-    a = -r['a1']; b = -r['a2']; c = -r['b2']
-    d = -r['a3']; k = -r['b3']; p = -r['c3']
-    sol1 = C1*exp(a*t)
-    sol2 = b*C1*exp(a*t)/(a-c) + C2*exp(c*t)
-    sol3 = C1*(d+b*k/(a-c))*exp(a*t)/(a-p) + k*C2*exp(c*t)/(c-p) + C3*exp(p*t)
-    return [Eq(x(t), sol1), Eq(y(t), sol2), Eq(z(t), sol3)]
 
 def _linear_3eq_order1_type2(x, y, z, t, r, eq):
     r"""

With the above patch test_ode still passes and this issue is solved:

In [1]: import sympy
   ...: from sympy import Eq, Derivative
   ...: a, b, c = sympy.symbols('a b c', cls=sympy.Function)
   ...: t, a_b, a_c = sympy.symbols('t a_b a_c', real=True)
   ...: 
   ...: eqs = [Eq(Derivative(a(t), t), (-a_b - a_c)*a(t)), Eq(Derivative(b(t), t), a_b*a(t)), Eq(Derivati
   ...: ve(c(t), t), a_c*a(t))]
   ...: sympy.dsolve(eqs)
   ...: 
   ...: 
Out[1]: 
⎡                        t⋅(-a_b - a_c)                              t⋅(-a_b - a_c)                    
⎢       -C₃⋅(a_b + a_c)⋅ℯ                                    C₃⋅a_b⋅ℯ                                t⋅
⎢a(t) = ────────────────────────────────, b(t) = C₁ + C₂⋅t + ──────────────────────, c(t) = C₂ + C₃⋅ℯ  
⎣                     a_c                                             a_c                              

            ⎤
(-a_b - a_c)⎥
            ⎥
            ⎦

I think the same goes for _linear_3eq_order1_type2 and _linear_3eq_order1_type3 (but not _linear_3eq_order1_typ4). These functions are all redundant in the face of _linear_neq_order1_type1.

@oscarbenjamin

This comment has been minimized.

Contributor

oscarbenjamin commented Nov 5, 2018

I just tried removing all 3 functions:

diff --git a/sympy/solvers/ode.py b/sympy/solvers/ode.py
index 8fc6f514a..f705f9ff5 100644
--- a/sympy/solvers/ode.py
+++ b/sympy/solvers/ode.py
@@ -7828,11 +7828,11 @@ def sysode_linear_3eq_order1(match_):
             if not j.has(x(t), y(t), z(t)):
                 raise NotImplementedError("Only homogeneous problems are supported, non-homogenous are not supported currently.")
     if match_['type_of_equation'] == 'type1':
-        sol = _linear_3eq_order1_type1(x, y, z, t, r, eq)
+        sol = _linear_neq_order1_type1(match_)
     if match_['type_of_equation'] == 'type2':
-        sol = _linear_3eq_order1_type2(x, y, z, t, r, eq)
+        sol = _linear_neq_order1_type1(match_)
     if match_['type_of_equation'] == 'type3':
-        sol = _linear_3eq_order1_type3(x, y, z, t, r, eq)
+        sol = _linear_neq_order1_type1(match_)
     if match_['type_of_equation'] == 'type4':
         sol = _linear_3eq_order1_type4(x, y, z, t, r, eq)
     if match_['type_of_equation'] == 'type6':

and all tests pass.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment