dsolve returns solution involving IOUSes (Integers of Unusual Size)

58 views
Skip to first unread message

G B

unread,
Apr 5, 2015, 12:27:07 AM4/5/15
to sy...@googlegroups.com
Still wrestling with dsolve...  Below is a call with an arbitrary differential equation.  Any idea why dsolve is returning terms with these enormous integers?  All of the coefficients are floats in this case.  The expression is impervious to .n() (as mentioned in my earlier question).  Converting to a numpy function to evaluate the results works, but throws an exception when called.

I can't seem to get past this point in the analysis.  Any idea how I can get this into a form I can continue working with?

A=symbols(r'A',cls=Function)
t=symbols(r't')
Eq4=-123456.78*A(t)-9876.54*A(t).diff(t)-0.00032*A(t).diff(t,2)+1357908.64
soln=dsolve(Eq4)
print(soln.n())

A(t) == C1*exp(125*t*(-33935533038108675 - sqrt(1151618536954453541512853661417481))/274877906944) + C2*exp(125*t*(-33935533038108675 + sqrt(1151618536954453541512853661417481))/274877906944) + 10.999060885923

fn=lambdify(t,soln.rhs,'numpy')

fn(3.2)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-42-bde0572cbbe1> in <module>()
----> 1 fn(3.2)

//anaconda/lib/python3.4/site-packages/numpy/__init__.py in <lambda>(_Dummy_73)

AttributeError: 'int' object has no attribute 'sqrt'

Ondřej Čertík

unread,
Apr 5, 2015, 12:40:49 AM4/5/15
to sympy
Hi,

I usually use symbols first, and only substitute into the final
result, that way I can see what's going on under the hood:

In [8]: var("a b c d")
Out[8]: (a, b, c, d)

In [9]: Eq4=a*A(t)-b*A(t).diff(t)-c*A(t).diff(t,2)+d

In [10]: Eq4
Out[10]:
2
d d
a⋅A(t) - b⋅──(A(t)) - c⋅───(A(t)) + d
dt 2
dt

In [11]: soln = dsolve(Eq4)

In [12]: soln
Out[12]:
⎛ ____________⎞ ⎛ ____________⎞
⎜ ╱ 2 ⎟ ⎜ ╱ 2 ⎟
t⋅⎝-b - ╲╱ 4⋅a⋅c + b ⎠ t⋅⎝-b + ╲╱ 4⋅a⋅c + b ⎠
──────────────────────── ────────────────────────
2⋅c 2⋅c d
A(t) = C₁⋅ℯ + C₂⋅ℯ - ─
a

In [13]: soln.rhs
Out[13]:
⎛ ____________⎞ ⎛ ____________⎞
⎜ ╱ 2 ⎟ ⎜ ╱ 2 ⎟
t⋅⎝-b - ╲╱ 4⋅a⋅c + b ⎠ t⋅⎝-b + ╲╱ 4⋅a⋅c + b ⎠
──────────────────────── ────────────────────────
2⋅c 2⋅c d
C₁⋅ℯ + C₂⋅ℯ - ─
a

In [14]: soln.rhs.subs({a: -123456.78, b: -9876.54, c: -0.00032, d:
1357908.64})Out[14]:
12.4999979732365⋅t -30864199.999998⋅t
C₁⋅ℯ + C₂⋅ℯ + 10.999060885923



Is the [14] what you wanted to get?

Ondrej
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/88a9432e-a322-460b-9556-ddd85cdad519%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

G B

unread,
Apr 5, 2015, 1:34:17 AM4/5/15
to sy...@googlegroups.com
Thanks, Ondrej.  In my case the coefficients are the result of a series of calculations earlier in the notebook.  They are symbolic up until the call to dsolve, but I am substituting values in right before the dsolve because it takes several minutes to return a solution otherwise (and gives a very complicated solution over a segmented domain).

I let it run to completion symbolically and tried substituting into the result and I do get something that looks more reasonable.  I can also run it through lambdify and plot the result.  So despite the slow computation, at least it looks like I can continue what I'm doing.

Any idea why the huge integers though?  I seem to remember somewhere that Sympy will try to convert floats to rationals under certain conditions.  Is that what's happening?  Can I turn that behavior off?

Any insight into why it won't evalf or lamdify if I substitute before dsolve rather than after?

Thanks, again!

Aaron Meurer

unread,
Apr 5, 2015, 2:39:52 AM4/5/15
to sy...@googlegroups.com
Looks like somewhere dsolve is converting the floats to rationals,
which end up being very large (large numerator / large denominator).
We should figure out where in dsolve this is happening and disable it.
It is likely coming from a call to solve(), since solve(rational=True)
does this.

As to it not evalfing, that looks like a separate bug.

Aaron Meurer
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/6f3572d3-07dd-40b3-a3a9-7673eddf1c87%40googlegroups.com.

Ondřej Čertík

unread,
Apr 5, 2015, 3:25:25 AM4/5/15
to sympy
On Sat, Apr 4, 2015 at 11:34 PM, G B <g.c.b....@gmail.com> wrote:
> Thanks, Ondrej. In my case the coefficients are the result of a series of
> calculations earlier in the notebook. They are symbolic up until the call
> to dsolve, but I am substituting values in right before the dsolve because
> it takes several minutes to return a solution otherwise (and gives a very
> complicated solution over a segmented domain).
>
> I let it run to completion symbolically and tried substituting into the
> result and I do get something that looks more reasonable. I can also run it
> through lambdify and plot the result. So despite the slow computation, at
> least it looks like I can continue what I'm doing.

Would you mind posting your whole calculation, if it is possible? We
are developing a very fast C++ core (https://github.com/sympy/csympy)
and I am always looking for examples of computations that people found
to be slow in SymPy, like in your case.

>
> Any idea why the huge integers though? I seem to remember somewhere that
> Sympy will try to convert floats to rationals under certain conditions. Is
> that what's happening? Can I turn that behavior off?
>
> Any insight into why it won't evalf or lamdify if I substitute before dsolve
> rather than after?

Aaron just answered these.

Ondrej

G B

unread,
Apr 5, 2015, 8:22:56 PM4/5/15
to sy...@googlegroups.com
Hey Ondrej--

I'll keep this in mind, but I don't think I can post these particular calculations as they're proprietary...

One thing this experience suggests though is that it might speed up dsolve if the coefficients are replaced with simple symbols such as you had and then resubstituted after the solve.  I don't know if that would be true in every case, or if it would prevent solutions in some cases...

Chris Smith

unread,
Apr 6, 2015, 9:16:47 AM4/6/15
to sy...@googlegroups.com
In line 4585 of ode.py RootOf instances return the Rational forms, e.g.,

>>> RootOf(2.3*x-2,0)
20/23. 

The default behavior of solve is "Float in, Float out" unless rational=True. (The docstring for solve should be updated to say that the default behavior is rational=None.)

Aaron Meurer

unread,
Apr 6, 2015, 12:13:34 PM4/6/15
to sy...@googlegroups.com
On Mon, Apr 6, 2015 at 8:16 AM, Chris Smith <smi...@gmail.com> wrote:
In line 4585 of ode.py RootOf instances return the Rational forms, e.g.,

>>> RootOf(2.3*x-2,0)
20/23. 

Ah. I suppose this will not be so easy to fix, as the float stuff with the polys does not work so well. We should at least not use RootOf unless we have to, though. 
 

The default behavior of solve is "Float in, Float out" unless rational=True. (The docstring for solve should be updated to say that the default behavior is rational=None.)

Good to know. I'm glad this is the behavior now. 

Aaron Meurer
 

On Saturday, April 4, 2015 at 11:27:07 PM UTC-5, G B wrote:
Still wrestling with dsolve...  Below is a call with an arbitrary differential equation.  Any idea why dsolve is returning terms with these enormous integers?  All of the coefficients are floats in this case.  The expression is impervious to .n() (as mentioned in my earlier question).  Converting to a numpy function to evaluate the results works, but throws an exception when called.

I can't seem to get past this point in the analysis.  Any idea how I can get this into a form I can continue working with?

A=symbols(r'A',cls=Function)
t=symbols(r't')
Eq4=-123456.78*A(t)-9876.54*A(t).diff(t)-0.00032*A(t).diff(t,2)+1357908.64
soln=dsolve(Eq4)
print(soln.n())

A(t) == C1*exp(125*t*(-33935533038108675 - sqrt(1151618536954453541512853661417481))/274877906944) + C2*exp(125*t*(-33935533038108675 + sqrt(1151618536954453541512853661417481))/274877906944) + 10.999060885923

fn=lambdify(t,soln.rhs,'numpy')

fn(3.2)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-42-bde0572cbbe1> in <module>()
----> 1 fn(3.2)

//anaconda/lib/python3.4/site-packages/numpy/__init__.py in <lambda>(_Dummy_73)

AttributeError: 'int' object has no attribute 'sqrt'

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.
Reply all
Reply to author
Forward
0 new messages