Skip to content
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.

Sign up for GitHub

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

Jump to bottom

Refactor liouville solver #21555

Merged
merged 4 commits into from Jun 3, 2021
Merged

Refactor liouville solver #21555

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
2d2bb9f
refactor liouville solver
Mohitbalwani26 Jun 1, 2021
0855ccf
remove unwanted variable
Mohitbalwani26 Jun 1, 2021
86d16ab
avoid repeated computation
Mohitbalwani26 Jun 1, 2021
72a2ecc
make dummy also class variable
Mohitbalwani26 Jun 1, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
3 changes: 2 additions & 1 deletion doc/src/modules/solvers/ode.rst
View file
Open in desktop
Expand Up @@ -111,7 +111,8 @@ Bernoulli

Liouville
^^^^^^^^^
.. autofunction:: sympy.solvers.ode.ode::ode_Liouville
.. autoclass:: sympy.solvers.ode.single::Liouville
:members:

Riccati_special_minus2
^^^^^^^^^^^^^^^^^^^^^^
Expand Down
94 changes: 2 additions & 92 deletions sympy/solvers/ode/ode.py
View file
Open in desktop
Expand Up @@ -991,7 +991,6 @@ class in it. Note that a hint may do this anyway if
a = Wild('a', exclude=[f(x)])
d = Wild('d', exclude=[df, f(x).diff(x, 2)])
e = Wild('e', exclude=[df])
k = Wild('k', exclude=[df])
n = Wild('n', exclude=[x, f(x), df])
c1 = Wild('c1', exclude=[x])
a3 = Wild('a3', exclude=[f(x), df, f(x).diff(x, 2)])
Expand Down Expand Up @@ -1055,6 +1054,7 @@ class in it. Note that a hint may do this anyway if
Factorable: ('factorable',),
RiccatiSpecial: ('Riccati_special_minus2',),
SecondNonlinearAutonomousConserved: ('2nd_nonlinear_autonomous_conserved',),
Liouville: ('Liouville',)
}
for solvercls in solvers:
solver = solvercls(ode)
Expand Down Expand Up @@ -1250,24 +1250,6 @@ def _powers(expr):


elif order == 2:
# Liouville ODE in the form
# f(x).diff(x, 2) + g(f(x))*(f(x).diff(x))**2 + h(x)*f(x).diff(x)
# See Goldstein and Braun, "Advanced Methods for the Solution of
# Differential Equations", pg. 98

s = d*f(x).diff(x, 2) + e*df**2 + k*df
r = reduced_eq.collect(f(x).diff(x)).match(s)
if r and r[d] != 0:
y = Dummy('y')
g = simplify(r[e]/r[d]).subs(f(x), y)
h = simplify(r[k]/r[d]).subs(f(x), y)
if y in h.free_symbols or x in g.free_symbols:
pass
else:
r = {'g': g, 'h': h, 'y': y}
matching_hints["Liouville"] = r
matching_hints["Liouville_Integral"] = r

# Homogeneous second order differential equation of the form
# a3*f(x).diff(x, 2) + b3*f(x).diff(x) + c3
# It has a definite power series solution at point x0 if, b3/a3 and c3/a3
Expand Down Expand Up @@ -3161,78 +3143,6 @@ def homogeneous_order(eq, *symbols):
return e


def ode_Liouville(eq, func, order, match):
r"""
Solves 2nd order Liouville differential equations.

The general form of a Liouville ODE is

.. math:: \frac{d^2 y}{dx^2} + g(y) \left(\!
\frac{dy}{dx}\!\right)^2 + h(x)
\frac{dy}{dx}\text{.}

The general solution is:

>>> from sympy import Function, dsolve, Eq, pprint, diff
>>> from sympy.abc import x
>>> f, g, h = map(Function, ['f', 'g', 'h'])
>>> genform = Eq(diff(f(x),x,x) + g(f(x))*diff(f(x),x)**2 +
... h(x)*diff(f(x),x), 0)
>>> pprint(genform)
2 2
/d \ d d
g(f(x))*|--(f(x))| + h(x)*--(f(x)) + ---(f(x)) = 0
\dx / dx 2
dx
>>> pprint(dsolve(genform, f(x), hint='Liouville_Integral'))
f(x)
/ /
| |
| / | /
| | | |
| - | h(x) dx | | g(y) dy
| | | |
| / | /
C1 + C2* | e dx + | e dy = 0
| |
/ /

Examples
========

>>> from sympy import Function, dsolve, Eq, pprint
>>> from sympy.abc import x
>>> f = Function('f')
>>> pprint(dsolve(diff(f(x), x, x) + diff(f(x), x)**2/f(x) +
... diff(f(x), x)/x, f(x), hint='Liouville'))
________________ ________________
[f(x) = -\/ C1 + C2*log(x) , f(x) = \/ C1 + C2*log(x) ]

References
==========

- Goldstein and Braun, "Advanced Methods for the Solution of Differential
Equations", pp. 98
- http://www.maplesoft.com/support/help/Maple/view.aspx?path=odeadvisor/Liouville

# indirect doctest

"""
# Liouville ODE:
# f(x).diff(x, 2) + g(f(x))*(f(x).diff(x, 2))**2 + h(x)*f(x).diff(x)
# See Goldstein and Braun, "Advanced Methods for the Solution of
# Differential Equations", pg. 98, as well as
# http://www.maplesoft.com/support/help/view.aspx?path=odeadvisor/Liouville
x = func.args[0]
f = func.func
r = match # f(x).diff(x, 2) + g*f(x).diff(x)**2 + h*f(x).diff(x)
y = r['y']
C1, C2 = get_numbered_constants(eq, num=2)
int = Integral(exp(Integral(r['g'], y)), (y, None, f(x)))
sol = Eq(int + C1*Integral(exp(-Integral(r['h'], x)), x) + C2, 0)
return sol


def ode_2nd_power_series_ordinary(eq, func, order, match):
r"""
Gives a power series solution to a second order homogeneous differential
Expand Down Expand Up @@ -6972,4 +6882,4 @@ def _nonlinear_3eq_order1_type5(x, y, z, t, eq):
#This import is written at the bottom to avoid circular imports.
from .single import (NthAlgebraic, Factorable, FirstLinear, AlmostLinear,
Bernoulli, SingleODEProblem, SingleODESolver, RiccatiSpecial,
SecondNonlinearAutonomousConserved, FirstExact)
SecondNonlinearAutonomousConserved, FirstExact, Liouville)
97 changes: 97 additions & 0 deletions sympy/solvers/ode/single.py
View file
Open in desktop
Expand Up @@ -980,5 +980,102 @@ def _get_general_solution(self, *, simplify: bool = True):
return [Eq(lhs, C2 + x), Eq(lhs, C2 - x)]


class Liouville(SinglePatternODESolver):
r"""
Solves 2nd order Liouville differential equations.

The general form of a Liouville ODE is

.. math:: \frac{d^2 y}{dx^2} + g(y) \left(\!
\frac{dy}{dx}\!\right)^2 + h(x)
\frac{dy}{dx}\text{.}

The general solution is:

>>> from sympy import Function, dsolve, Eq, pprint, diff
>>> from sympy.abc import x
>>> f, g, h = map(Function, ['f', 'g', 'h'])
>>> genform = Eq(diff(f(x),x,x) + g(f(x))*diff(f(x),x)**2 +
... h(x)*diff(f(x),x), 0)
>>> pprint(genform)
2 2
/d \ d d
g(f(x))*|--(f(x))| + h(x)*--(f(x)) + ---(f(x)) = 0
\dx / dx 2
dx
>>> pprint(dsolve(genform, f(x), hint='Liouville_Integral'))
f(x)
/ /
| |
| / | /
| | | |
| - | h(x) dx | | g(y) dy
| | | |
| / | /
C1 + C2* | e dx + | e dy = 0
| |
/ /

Examples
========

>>> from sympy import Function, dsolve, Eq, pprint
>>> from sympy.abc import x
>>> f = Function('f')
>>> pprint(dsolve(diff(f(x), x, x) + diff(f(x), x)**2/f(x) +
... diff(f(x), x)/x, f(x), hint='Liouville'))
________________ ________________
[f(x) = -\/ C1 + C2*log(x) , f(x) = \/ C1 + C2*log(x) ]

References
==========

- Goldstein and Braun, "Advanced Methods for the Solution of Differential
Equations", pp. 98
- http://www.maplesoft.com/support/help/Maple/view.aspx?path=odeadvisor/Liouville

# indirect doctest

"""
hint = "Liouville"
has_integral = True
order = [2]

def _wilds(self, f, x, order):
d = Wild('d', exclude=[f(x).diff(x), f(x).diff(x, 2)])
e = Wild('e', exclude=[f(x).diff(x)])
k = Wild('k', exclude=[f(x).diff(x)])
return d, e, k

def _equation(self, fx, x, order):
# Liouville ODE in the form
# f(x).diff(x, 2) + g(f(x))*(f(x).diff(x))**2 + h(x)*f(x).diff(x)
# See Goldstein and Braun, "Advanced Methods for the Solution of
# Differential Equations", pg. 98
d, e, k = self.wilds()
return d*fx.diff(x, 2) + e*fx.diff(x)**2 + k*fx.diff(x)

def _verify(self, fx):
d, e, k = self.wilds_match()
y = Dummy('y')
x = self.ode_problem.sym
self.g = simplify(e/d).subs(fx, y)
self.h = simplify(k/d).subs(fx, y)
if y in self.h.free_symbols or x in self.g.free_symbols:
return False
return True

def _get_general_solution(self, *, simplify: bool = True):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What does simplify: bool=True mean? And how is it used (whatever the name) in this function? Perhaps it should just be dropped.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess this was used to simplify the solution, not sure though. @oscarbenjamin can you please explain a little?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't see the simplify flag being used anywhere. Not sure what it means. Maybe it should be dropped. In dsolve simplification is controlled elsewhere...

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How do you read that syntax? what is 'simplify: bool=True'?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, static type-check reads as "simplify is a bool and its value is True by default". Let's drop this arg -- there is no simplification being requested or done by any of these functions.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It might be needed for some solvers so I suggest to finish the refactoring before deciding about that. These classes and methods are all internal so it should be safe to delete the argument later.

d, e, k = self.wilds_match()
fx = self.ode_problem.func
x = self.ode_problem.sym
y = Dummy('y')
C1, C2 = self.ode_problem.get_numbered_constants(num=2)
int = Integral(exp(Integral(self.g, y)), (y, None, fx))
gen_sol = Eq(int + C1*Integral(exp(-Integral(self.h, x)), x) + C2, 0)

return [gen_sol]


# Avoid circular import:
from .ode import dsolve