Add meijerg = False in dsolve

43 views
Skip to first unread message

Tom van Woudenberg

unread,
Jan 17, 2023, 9:57:02 AM1/17/23
to sympy
Hi there,

When solving differential equations with Heaviside functions, a meijerg function appears in the outcome, while a (simpler) expression is possible as well. For example a beam problem:

import sympy as sp
w = sp.symbols('w', cls=sp.Function)
x = sp.symbols('x')
C1, C2, C3, C4 = sp.symbols('C1 C2 C3 C4')
Q, a, b, EI = sp.symbols('Q, a, b, EI')
q = Q*(1 - sp.Heaviside(x-a))+F*sp.DiracDelta(x-(a+b))
diffeq = sp.Eq(EI*sp.diff(w(x),x,4),q)
w = sp.dsolve(diffeq)
w = w.rhs
print(w)

gives:
C1 + x**3*(C4 - Q*Piecewise((Piecewise((0, Abs(x/a) < 1), (a*meijerg(((2, 1), ()), ((), (1, 0)), x/a), True)), a > 0), (Integral(Heaviside(-a + x), x), True))/(6*EI) + 35*Heaviside(-a - b + x)/(6*EI)) + x**2*(C3 + Q*Piecewise((Piecewise((0, Abs(x/a) < 1), (a**2*meijerg(((3, 1), ()), ((), (2, 0)), x/a), True)), a > 0), (Integral(x*Heaviside(-a + x), x), True))/(2*EI) - 35*a*Heaviside(-a - b + x)/(2*EI) - 35*b*Heaviside(-a - b + x)/(2*EI)) + x*(C2 - Q*Piecewise((Piecewise((0, Abs(x/a) < 1), (a**3*meijerg(((4, 1), ()), ((), (3, 0)), x/a), True)), a > 0), (Integral(x**2*Heaviside(-a + x), x), True))/(2*EI) + 35*(a + b)**2*Heaviside(-a - b + x)/(2*EI)) + Q*x**4/(24*EI) + Q*Piecewise((Piecewise((0, Abs(x/a) < 1), (a**4*meijerg(((5, 1), ()), ((), (4, 0)), x/a), True)), a > 0), (Integral(x**3*Heaviside(-a + x), x), True))/(6*EI) - 35*(a + b)**3*Heaviside(-a - b + x)/(6*EI)
Thus including meijerg functions.

If the same result is found by integrating multiple times with meijerg=False:
import sympy as sp
w = sp.symbols('w', cls=sp.Function)
x = sp.symbols('x')
C1, C2, C3, C4 = sp.symbols('C1 C2 C3 C4')
Q, a, b, EI = sp.symbols('Q, a, b, EI')
q = Q*(1 - sp.Heaviside(x-a))+F*sp.DiracDelta(x-(a+b))
V = sp.integrate(-q,x,meijerg=False)+C1
M = sp.integrate(V,x,meijerg=False)+C2
kappa = M / EI
phi = sp.integrate(kappa,x,meijerg=False)+C3
w = sp.integrate(-phi,x,meijerg=False)+C4
print(w)


the following result is obtained:
-C3*x + C4 - (C1*x**3/6 + C2*x**2/2 - Q*(x**4/24 - (a**4/24 - a**3*x/6 + a**2*x**2/4 - a*x**3/6 + x**4/24)*Heaviside(-a + x)) - 35*(-a*x**2/2 + a*x*(a + b) - a*(a + b)**2/2 - b*x**2/2 + b*x*(a + b) - b*(a + b)**2/2 + x**3/6 - x*(a + b)**2/2 + (a + b)**3/3)*Heaviside(-a - b + x))/EI

Would it be possible to implement a meijerg=false option for dsolve as well?

Kind regards,
Tom van Woudenberg

Oscar Benjamin

unread,
Jan 17, 2023, 1:17:02 PM1/17/23
to sy...@googlegroups.com
On Tue, 17 Jan 2023 at 14:57, 'Tom van Woudenberg' via sympy
<sy...@googlegroups.com> wrote:
>
> Would it be possible to implement a meijerg=false option for dsolve as well?

It would be better to have a simple way of disabling integration in dsolve.

Jeremy Monat

unread,
Jan 17, 2023, 1:22:24 PM1/17/23
to sy...@googlegroups.com

I recall Aaron Meurer mentioned that it would be nice to have a simpler way of setting that flag.

Jeremy Monat

--
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 view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxS99wYz%2B0Kw%2BByW9Y0M%3DC-v0Y%2Bj%3DiDZ9KcWe-CU6s1Lww%40mail.gmail.com.

Aaron Meurer

unread,
Jan 17, 2023, 1:36:19 PM1/17/23
to sy...@googlegroups.com
Yes. I was going to mention that to do this, you should disable integration. Right now, the only way to do this is to get the specific hint and add "_Integral" to the end. Unfortunately, classify_ode hangs on this differential equation (it gets stuck in the solve() call in the nth algebraic matcher), so it's hard to say which hint to use. I agree that integrate=False ought to just be a flag to dsolve(). This shouldn't be a difficult change to make.

If you were to do this, you could just call doit(meijerg=False) on the resulting expression.

Aaron Meurer

Jeremy Monat

unread,
Jan 17, 2023, 1:38:48 PM1/17/23
to sy...@googlegroups.com
The simplest way to disable integration is with the all_Integral hint because you do not need to know which hint to supply: for any hint with a corresponding _Integral hint, all_Integral only returns the _Integral hint.

Jeremy Monat

Tom van Woudenberg

unread,
Jan 18, 2023, 3:33:28 AM1/18/23
to sympy
Disabling integration would be a neat solution! However, as Aaron stated, dsolve() hangs when using all_Integral. What would be the best way to proceed? Would a flag integrate=False be an appropriate issue to start?

Op dinsdag 17 januari 2023 om 19:38:48 UTC+1 schreef Jeremy Monat:

Oscar Benjamin

unread,
Jan 18, 2023, 6:47:05 AM1/18/23
to sy...@googlegroups.com
On Wed, 18 Jan 2023 at 08:33, 'Tom van Woudenberg' via sympy
<sy...@googlegroups.com> wrote:
>
> Disabling integration would be a neat solution! However, as Aaron stated, dsolve() hangs when using all_Integral.

That's a separate issue that should also be fixed.

> What would be the best way to proceed? Would a flag integrate=False be an appropriate issue to start?

Yes.

--
Oscar

Tom van Woudenberg

unread,
Jan 18, 2023, 9:41:56 AM1/18/23
to sympy
I opened two issues: https://github.com/sympy/sympy/issues/24546 and https://github.com/sympy/sympy/issues/24547

Op woensdag 18 januari 2023 om 12:47:05 UTC+1 schreef Oscar:
Reply all
Reply to author
Forward
0 new messages