Problems with dsolve
Francesco Paparella
I'm pretty new to sympy, which I like a lot: it feels like a very
pythonic approach to symbolic computations.
However, I've just had a rough encounter with 'dsolve' this afternoon.
There's a mishandling of 1st-order constant coefficient equations.
Is it a known problem? (I'm using v.0.6.7 on ubuntu).
A sample session follows.
Francesco
-----------------------------------------------------------------------
Python 2.6.5 (r265:79063, Apr 16 2010, 13:09:56)
[GCC 4.4.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> from sympy import *
>>> a=Symbol('a')
>>> b=Symbol('b')
>>> t=Symbol('t')
>>> f=Function('f')
>>>
>>> dsolve(f(t).diff(t) - f(t), f(t))
f(t) == exp(C1 + t)
>>> #right... but what if I have an initial condition f(0)<0 ?
...
>>> dsolve(f(t).diff(t) - a - f(t), f(t))
f(t) == -a - exp(C1 + t)
>>> #so far, so good (except for f(0)>-a)
...
>>> dsolve(f(t).diff(t) - a - b*f(t), f(t))
f(t) == -(a + exp(C1 + b*t))/b
>>> #just like above
...
>>> dsolve(f(t).diff(t) - a - (b-a)*f(t), f(t))
f(t) == (1 + a*exp(C1 + a*t - b*t))/(a*exp(C2 + a*t - b*t) - b*exp(C3 +
a*t - b*t))
>>> #what? Three constants of integration?
...
>>> dsolve(f(t).diff(t) - a - (b-a)*f(t), f(t), hint='1st_linear')
f(t) == exp(-t*(a - b))*(C1 + Integral(a*exp(t*(a - b)), t))
>>> #Hey! this is correct (I can even impose any initial condition), but
why it doesn't carry out the integral?
...
>>> simplify(dsolve(f(t).diff(t) - a - (b-a)*f(t), f(t),
hint='1st_linear').rhs)
exp(b*t - a*t)*Integral(a*exp(a*t)*exp(-b*t), t) + C1*exp(b*t - a*t)
>>> #do the damn integral, would you?
...
>>> integrate(a*exp(t*(a - b)), t)
Integral(a*exp(t*(a - b)), t)
>>> #oh... so, maybe you just can't integrate?
...
>>> integrate(a*exp(t*b), t)
a*exp(b*t)/b
>>> #yes, you can! What's the problem? You don't like parentheses?
...
>>> import sympy
>>> sympy.__version__
'0.6.7'
smichr
> There's a mishandling of 1st-order constant coefficient equations.
> Is it a known problem? (I'm using v.0.6.7 on ubuntu).
> A sample session follows.
ring in on this, but regarding sympy not doing the integral...you can
help it out by simply using expand on your function:
>>> integrate(a*exp(t*(b-a)).expand(), t)
a*exp(b*t)/(b*exp(a*t) - a*exp(a*t))
/c
smichr
>
> f(t) == (1 + a*exp(C1 + a*t - b*t))/(a*exp(C2 + a*t - b*t) - b*exp(C3 +
> a*t - b*t))>>> #what? Three constants of integration?
h[8] >>> dsolve(f(t).diff(t) - a - (b-a)*f(t), f(t))
b*exp(C3 + a*t - b*
t))
h[8] >>> n/powsimp(factor(d))
(1 + a*exp(C1 + a*t - b*t))*exp(b*t - a*t)/(a*exp(C2) - b*exp(C3))
h[9] >>> _.subs(1/fraction(_)[1],Symbol('C4'))
C4*(1 + a*exp(C1 + a*t - b*t))*exp(b*t - a*t)
/c
Franz
That's very useful. I guess 'expand everything' is the way to go in
sympy, isn't it?
> The integration is undergoing a big overhaul and Aaron will probably
> ring in on this, but regarding sympy not doing the integral...you can
> help it out by simply using expand on your function:
>
> >>> integrate(a*exp(t*(b-a)).expand(), t)
> a*exp(b*t)/(b*exp(a*t) - a*exp(a*t))
r=a*exp(b*t)/(b*exp(a*t) - a*exp(a*t))
what is the sympytonic way to go back to an expression with a single
exponential?
I've found:
powsimp(collect(r, exp(a*t)))
By the way, 'collect' doesn't seem to work inside exponentials:
collect( powsimp(collect(r, exp(a*t))), t)
does not collect the 't'.
Franz
By the way, when I attempt:
factor(d)
I get:
PolynomialError: Can't decompose exp(C3)
Is there something wrong with my installation, or am I using an old
version (0.6.7)?
Aaron S. Meurer
On Feb 2, 2011, at 12:00 PM, Francesco Paparella wrote:
> Hello,
>
> I'm pretty new to sympy, which I like a lot: it feels like a very
> pythonic approach to symbolic computations.
>
> However, I've just had a rough encounter with 'dsolve' this afternoon.
> There's a mishandling of 1st-order constant coefficient equations.
> Is it a known problem? (I'm using v.0.6.7 on ubuntu).
> A sample session follows.
>
> Francesco
>
> -----------------------------------------------------------------------
> Python 2.6.5 (r265:79063, Apr 16 2010, 13:09:56)
> [GCC 4.4.3] on linux2
> Type "help", "copyright", "credits" or "license" for more information.
>>>> from sympy import *
>>>> a=Symbol('a')
>>>> b=Symbol('b')
>>>> t=Symbol('t')
>>>> f=Function('f')
>>>>
>>>> dsolve(f(t).diff(t) - f(t), f(t))
> f(t) == exp(C1 + t)
>>>> #right... but what if I have an initial condition f(0)<0 ?
For now, you will have to solve for the constants manually if you have an initial condition, because they aren't implanted yet. See
> ...
>>>> dsolve(f(t).diff(t) - a - f(t), f(t))
> f(t) == -a - exp(C1 + t)
>>>> #so far, so good (except for f(0)>-a)
> ...
>>>> dsolve(f(t).diff(t) - a - b*f(t), f(t))
> f(t) == -(a + exp(C1 + b*t))/b
>>>> #just like above
Wait, what is the problem here?
> ...
>>>> dsolve(f(t).diff(t) - a - (b-a)*f(t), f(t))
> f(t) == (1 + a*exp(C1 + a*t - b*t))/(a*exp(C2 + a*t - b*t) - b*exp(C3 +
> a*t - b*t))
>>>> #what? Three constants of integration?
So here is what happens. dsolve has an internal function called constantsimp that simplifies constants of integration in the solution. So instead of getting something like f(x) == 2*C1*x, you just get f(x) == C1*x. The problem is that sometimes after solving for f(x), the constants split up, and the constantsimp simplifies (and renumbers) them into two constants. The solution is still correct except that the constants will not be independent of each other. From the docstring of dsolve():
``simplify`` enables simplification by odesimp(). See its
docstring for more information. Turn this off, for example,
to disable solving of solutions for func or simplification
of arbitrary constants. It will still integrate with this
hint. Note that the solution may contain more arbitrary
constants than the order of the ODE with this option
enabled.
So the solution is to add simplify=False to dsolve(), which will disable all simplification, including constant simplification and solving for f(x) (you will have to do those manually). Doing this gives:
In [12]: dsolve(f(t).diff(t) - a - (b-a)*f(t), f(t), simplify=False)
Out[12]:
log(-a + (a - b)⋅f(t))
────────────────────── = C₁ - t
a - b
Notice how when you solve this for f(t), it splits C1 into multiple places:
In [14]: solve(dsolve(f(t).diff(t) - a - (b-a)*f(t), f(t), simplify=False), f(t))
Out[14]:
⎡ C₁⋅b + a⋅t - C₁⋅a - b⋅t ⎤
⎢ 1 + a⋅ℯ ⎥
⎢───────────────────────────────────────────────────────⎥
⎢ C₁⋅b + a⋅t - C₁⋅a - b⋅t C₁⋅b + a⋅t - C₁⋅a - b⋅t⎥
⎣a⋅ℯ - b⋅ℯ ⎦
constantsimp() takes a cautious approach and renumbers each of these C1 differently, because for example C1*b != C1*a.
> ...
>>>> dsolve(f(t).diff(t) - a - (b-a)*f(t), f(t), hint='1st_linear')
> f(t) == exp(-t*(a - b))*(C1 + Integral(a*exp(t*(a - b)), t))
>>>> #Hey! this is correct (I can even impose any initial condition), but
> why it doesn't carry out the integral?
> ...
>>>> simplify(dsolve(f(t).diff(t) - a - (b-a)*f(t), f(t),
> hint='1st_linear').rhs)
> exp(b*t - a*t)*Integral(a*exp(a*t)*exp(-b*t), t) + C1*exp(b*t - a*t)
>>>> #do the damn integral, would you?
> ...
>>>> integrate(a*exp(t*(a - b)), t)
> Integral(a*exp(t*(a - b)), t)
>>>> #oh... so, maybe you just can't integrate?
> ...
>>>> integrate(a*exp(t*b), t)
> a*exp(b*t)/b
>>>> #yes, you can! What's the problem? You don't like parentheses?
Yes, so you are discovering one of the (presently) weaker points of sympy. Chris's workarounds should work for now. In the future, the Risch Algorithm will be able to handle this just fine. For example, from my integration3 branch:
In [15]: risch_integrate(a*exp(t*(a - b)), t)
Out[15]:
t⋅(a - b)
a⋅ℯ
────────────
a - b
but the other version points out a bug in risch_integrate(), which I will look into immediately.
> ...
>>>> import sympy
>>>> sympy.__version__
> '0.6.7'
Aaron Meurer
Aaron S. Meurer
Factor in the current version is very limited (only actual polynomials allowed), but this is all fixed in the current development branch and will be fixed in the current master.
Aaron Meurer
Aaron S. Meurer
powsimp is the function that should recombine exponents. You will need to factor the expression so that they are actually multiplying each other first. In the current development branch---and in the next release---this works:
In [22]: powsimp(factor(r))
Out[22]:
b⋅t - a⋅t
-a⋅ℯ
─────────────
a - b
By the way, if you don't want to wait for the next release, it is very easy to work off of the development branch. See https://github.com/sympy/sympy/wiki/Getting-The-Bleeding-Edge for a little guide.
Aaron Meurer
Franz
first of all, thanks for such an extensive and informative answer.
> >>>> dsolve(f(t).diff(t) - f(t), f(t))
> > f(t) == exp(C1 + t)
> >>>> #right... but what if I have an initial condition f(0)<0 ?
>
> For now, you will have to solve for the constants manually if you have an initial condition, because they aren't implanted yet.
that dsolve has only found some solutions, not the general integral.
In fact:
f(0)=exp(C1)
but the exponential is always positive, therefore the solutions with
f(0)<=0 are not representable with the formula given by dsolve.
Of course, at calculus 101 we all defined C2=exp(C1), and then let C2
have any real value. But formally this is wrong (exponentials are
positive). Another way to say the same thing, is that the integral of
df/f is log(|f|), while the default solver (separation of variables?)
seems to use df/f = d log(f) (which is not ok for f<=0).
In fact, the solver does the right thing, when given the hint
'1st_linear':
>>> dsolve(f(t).diff(t) - f(t), f(t), hint='1st_linear')
f(t) == C1*exp(t)
You will agree this is quite more general than f(t) == exp(C1 + t)
> So here is what happens. dsolve has an internal function called constantsimp that simplifies constants of integration in the solution. So instead of getting something like f(x) == 2*C1*x, you just get f(x) == C1*x. The problem is that sometimes after solving for f(x), the constants split up, and the constantsimp simplifies (and renumbers) them into two constants. The solution is still correct except that the constants will not be independent of each other. From the docstring of dsolve():
>
> ``simplify`` enables simplification by odesimp(). See its
> docstring for more information. Turn this off, for example,
> to disable solving of solutions for func or simplification
> of arbitrary constants. It will still integrate with this
> hint. Note that the solution may contain more arbitrary
> constants than the order of the ODE with this option
> enabled.
>
> So the solution is to add simplify=False to dsolve(), which will disable all simplification, including constant simplification and solving for f(x) (you will have to do those manually). Doing this gives:
>
> In [12]: dsolve(f(t).diff(t) - a - (b-a)*f(t), f(t), simplify=False)
> Out[12]:
> log(-a + (a - b)⋅f(t))
> ────────────────────── = C₁ - t
> a - b
>
> Notice how when you solve this for f(t), it splits C1 into multiple places:
>
> In [14]: solve(dsolve(f(t).diff(t) - a - (b-a)*f(t), f(t), simplify=False), f(t))
> Out[14]:
> ⎡ C₁⋅b + a⋅t - C₁⋅a - b⋅t ⎤
> ⎢ 1 + a⋅ℯ ⎥
> ⎢───────────────────────────────────────────────────────⎥
> ⎢ C₁⋅b + a⋅t - C₁⋅a - b⋅t C₁⋅b + a⋅t - C₁⋅a - b⋅t⎥
> ⎣a⋅ℯ - b⋅ℯ ⎦
>
> constantsimp() takes a cautious approach and renumbers each of these C1 differently, because for example C1*b != C1*a.
this 'cautious approach' is non very cautious and that it should be
reconsidered in the future? Would it be so terrible to never absorb
constants into an integration constant?
> Yes, so you are discovering one of the (presently) weaker points of sympy. Chris's workarounds should work for now. In the future, the Risch Algorithm will be able to handle this just fine. For example, from my integration3 branch:
>
> In [15]: risch_integrate(a*exp(t*(a - b)), t)
> Out[15]:
> t⋅(a - b)
> a⋅ℯ
> ────────────
> a - b
F
Aaron S. Meurer
On Feb 2, 2011, at 4:50 PM, Franz wrote:
> Hi Aaron,
>
> first of all, thanks for such an extensive and informative answer.
Well, that's because I wrote the ode module :)
>
>>>>>> dsolve(f(t).diff(t) - f(t), f(t))
>>> f(t) == exp(C1 + t)
>>>>>> #right... but what if I have an initial condition f(0)<0 ?
>>
>> For now, you will have to solve for the constants manually if you have an initial condition, because they aren't implanted yet.
>
> I'm not complaining about manually solving for C1. I'm trying to say
> that dsolve has only found some solutions, not the general integral.
> In fact:
> f(0)=exp(C1)
> but the exponential is always positive, therefore the solutions with
> f(0)<=0 are not representable with the formula given by dsolve.
> Of course, at calculus 101 we all defined C2=exp(C1), and then let C2
> have any real value. But formally this is wrong (exponentials are
> positive). Another way to say the same thing, is that the integral of
> df/f is log(|f|), while the default solver (separation of variables?)
> seems to use df/f = d log(f) (which is not ok for f<=0).
> In fact, the solver does the right thing, when given the hint
> '1st_linear':
>>>> dsolve(f(t).diff(t) - f(t), f(t), hint='1st_linear')
> f(t) == C1*exp(t)
> You will agree this is quite more general than f(t) == exp(C1 + t)
Yes, I was going to suggest that you use the hints system. The way that the solver works is that it just recognizes a form and then just plugs the coefficients into a solution, so the solution given by one might not be as pleasing in one respect or another as the solution given by another. Usually when I use dsolve() I use the 'all' hint and then look through to see which solution looks the best to me (the 'best' hint usually isn't too bad at picking it out either, actually). It's slower because it solves the ode using every matching method, but if you are only solving one or two odes, it isn't bad and worth it.
Actually, imo we should fix constantsimp to reduce that solution to C1*exp(t), because that is a much better looking solution and the separable one is what is given by default.
Well, if you don't like it, you can always just pass the simplify=False flag to dsolve(). Perhaps we could make it more cautious and have it not simplify the constants if the result produces more constants than the order of the ode, or maybe we could only constantsimp() before solving for f(x) (though this last option will miss even some cases where it won't produce more than one constant).
I wrote the constantsimp routine because I was annoyed how Maple would return things like
> dsolve(diff(f(x), x, x) + f(x) - 2 + sin(x));
1 1
f(x) = sin(x) _C2 + cos(x) _C1 + 2 - - sin(x) + - cos(x) x
2 2
where the -sin(x)/2 term is completely redundant because of the sin(x)*_C2 term. SymPy would do the same thing (this is an artifact of the variation of parameters method), if it weren't for constantsimp(), and also the undetermined coefficients method that SymPy implements and Maple does not, which completely avoids this.
>
>> Yes, so you are discovering one of the (presently) weaker points of sympy. Chris's workarounds should work for now. In the future, the Risch Algorithm will be able to handle this just fine. For example, from my integration3 branch:
>>
>> In [15]: risch_integrate(a*exp(t*(a - b)), t)
>> Out[15]:
>> t⋅(a - b)
>> a⋅ℯ
>> ────────────
>> a - b
>
> Wonderful! Sympy is on the road for world domination!
> F
Well, this isn't a particularly tricky integral, but I agree :)
Aaron Meurer
Chris Smith
>
> factor(d)
>
> I get:
>
> PolynomialError: Can't decompose exp(C3)
I was using my casy branch at github (smichr acct).
/c
Chris Smith
>
> factor(d)
>
> I get:
>
> PolynomialError: Can't decompose exp(C3)
I was using my casy branch at github (smichr acct).
/c
Vinzent Steinberg
> That's very useful. I guess 'expand everything' is the way to go in
> sympy, isn't it?
and sometimes it makes your expressions so huge that you can't
simplify them properly afterwards.
Vinzent