Bernoulli ODE has one case where it should be solved using different
algorithm. Currently we don't detect it and hilarity ensues in the form
of divide by zero :)

https://code.plaes.org/git/sympy.git/commit/?h=patches-0514&id=56524bf458074ac9635f537ef6f83c4727ea7a56



--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups 
"sympy-patches" group.
To post to this group, send email to sympy-patches@googlegroups.com
To unsubscribe from this group, send email to 
sympy-patches+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sympy-patches?hl=en
-~----------~----~----~----~------~----~------~--~---

>From ffc4f78386316e5a876d7d83672b72319fa01140 Mon Sep 17 00:00:00 2001
From: Priit Laes <pl...@plaes.org>
Date: Sat, 11 Apr 2009 21:15:13 +0300
Subject: [PATCH 2/2] Added corner case for Bernoulli equation.

Signed-off-by: Priit Laes <pl...@plaes.org>
---
 sympy/solvers/solvers.py            |    8 ++++++--
 sympy/solvers/tests/test_solvers.py |    1 +
 2 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index 16a14c2..e6928dc 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -569,6 +569,7 @@ def solve_ODE_first_order(eq, f):
     and Bernoulli cases are implemented.
     """
     from sympy.integrals.integrals import integrate
+    C1 = Symbol("C1")
     x = f.args[0]
     f = f.func
 
@@ -581,16 +582,19 @@ def solve_ODE_first_order(eq, f):
     if r:
         t = exp(integrate(r[b]/r[a], x))
         tt = integrate(t*(-r[c]/r[a]), x)
-        return (tt + Symbol("C1"))/t
+        return (tt + C1)/t
 
     #Bernoulli case: a(x)*f'(x)+b(x)*f(x)+c(x)*f(x)^n = 0
     n = Wild('n', exclude=[f(x)])
 
     r = eq.match(a*diff(f(x),x) + b*f(x) + c*f(x)**n)
     if r:
+        if r[n] == 1:
+            return C1*exp(integrate(-(r[b]+r[c]), x))
+        # r[n] != 1 ie, the real bernoulli case
         t = exp((1-r[n])*integrate(r[b]/r[a],x))
         tt = (r[n]-1)*integrate(t*r[c]/r[a],x)
-        return ((tt + Symbol("C1"))/t)**(1/(1-r[n]))
+        return ((tt + C1)/t)**(1/(1-r[n]))
 
     #other cases of first order odes will be implemented here
 
diff --git a/sympy/solvers/tests/test_solvers.py b/sympy/solvers/tests/test_solvers.py
index 8d02fb7..a4d9d1f 100644
--- a/sympy/solvers/tests/test_solvers.py
+++ b/sympy/solvers/tests/test_solvers.py
@@ -158,6 +158,7 @@ def test_ODE_first_order():
     assert dsolve(3*f(x).diff(x) -1, f(x)) == x/3 + C1
     assert dsolve(x*f(x).diff(x) -1, f(x)) == log(x) + C1
     assert dsolve(x*f(x).diff(x)+f(x)-f(x)**2,f(x)) == 1/(x*(C1 + 1/x))
+    assert dsolve(f(x).diff(x)+x*f(x)-f(x),f(x)) == C1*exp(x - x**2/2)
 
 def test_ODE_second_order():
     f = Function('f')
-- 
1.6.3

Reply via email to