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

[GSoC] Kovacic Solver - Part 1 #21770

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

[GSoC] Kovacic Solver - Part 1 #21770

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
67 changes: 67 additions & 0 deletions sympy/solvers/ode/kovacic.py
View file
Open in desktop
@@ -0,0 +1,67 @@
from .riccati import solve_riccati

from sympy import S, Eq, exp, Integral, Derivative, Function
from sympy.core.symbol import symbols, Wild
from sympy.solvers.ode.ode import get_numbered_constants, _remove_redundant_solutions, constantsimp
from sympy.solvers.ode.subscheck import checkodesol

def match_2nd_order(eq, z, x):
a = Wild('a', exclude=[z, z.diff(x), z.diff(x)*2])
b = Wild('b', exclude=[z, z.diff(x), z.diff(x)*2])
c = Wild('c', exclude=[z, z.diff(x), z.diff(x)*2])
match = eq.match(a*z.diff(x, 2) + b*z.diff(x) + c*z)

if a not in match or b not in match or c not in match:
raise ValueError("Invalid Second Order Linear Homogeneous Equation")
return match[a], match[b], match[c]
Copy link
Contributor

Choose a reason for hiding this comment

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

Shouldn't this check that the matches are rational functions?


def find_kovacic_simple(x, a, b, c):
# Find solution for y(x).diff(x, 2) = r(x)*y(x)
r = (b**2 + 2*a*b.diff(x) - 2*a.diff(x)*b - 4*a*c)/4*a**2

# Riccati equation to be solved is g(x).diff(x) + g(x)**2 - r
# Comparing with f(x).diff(x) = b_0(x) + b_1(x)*f(x) + b_2(x)*f(x)**2
# we can see that b_0 = r, b_1 = 0, b_2 = -1
b0, b1, b2 = r, S(0), -S(1)
g = Function('g')
ric_sol = [sol.rhs for sol in solve_riccati(g(x), x, b0, b1, b2)]

C1 = symbols('C1')
return set(map(lambda sol: exp(Integral(sol.subs(C1, 0), x).doit()), ric_sol))
Copy link
Contributor

Choose a reason for hiding this comment

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

This should use Dummy so it doesn't clash with symbols from other calls to dsolve

Copy link
Contributor

Choose a reason for hiding this comment

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

Also don't use map/lambda like this. A comprehension is clearer:

return {exp(Integral(sol.subs(C1, 0), x).doit()) for sol in ric_sol}

Copy link
Member

Choose a reason for hiding this comment

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

Do we have the same problem with Function('g') above? Will this code break if the ODE already has a function named g in it?

Copy link
Member

Choose a reason for hiding this comment

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

If I read this correctly, it is assuming that C1 is the constant that comes from the recursive dsolve solution. Is this something that is always true? How do we know the constant is 0? Can this be expressed as an initial condition?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

If I read this correctly, it is assuming that C1 is the constant that comes from the recursive dsolve solution. Is this something that is always true? How do we know the constant is 0? Can this be expressed as an initial condition?

Sometimes the riccati solver returns a general solution. However, we only need one particular solution for the auxiliary riccati equation to find the general solution for the second order ODE. That's the reason I'm substituting it with 0. We could substitute it with anything finite, it wouldn't matter.

Copy link
Contributor

Choose a reason for hiding this comment

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

If we only need a particular solution is there a way of calling the Ricatti solver that just tried to find the first particular solution as quickly as possible?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Actually the code is supposed to only find particular solutions. However what happens sometimes is that some coefficients remain undetermined when the undetermined coefficients are solved for the auxiliary equation, and so a general solution may appear instead of a particular solution. In either case, the time to calculate the solution doesn't change and it doesn't involve any slow steps like calculating any integrals.

Copy link
Contributor

Choose a reason for hiding this comment

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

Does it just return the first particular solution found or does it still run through finding all of them?


def find_kovacic_sol(eq):
z = list(eq.atoms(Derivative))[0].args[0]
x = list(z.free_symbols)[0]
eq = eq.expand().collect(z)
a, b, c = match_2nd_order(eq, z, x)

# Transform the differential equation to a simpler form
# using z(x) = y(x)*exp(Integral(-b/(2*a))) and find its solution
ysol = find_kovacic_simple(x, a, b, c)

zsol = list(map(lambda sol: sol*exp(Integral(-b/(2*a), x).doit()), ysol))
zsol = _remove_redundant_solutions(Eq(eq, 0), list(map(lambda sol: Eq(z, sol), zsol)), 2, x)

C1, C2 = symbols('C1 C2')
print(zsol)
if len(zsol) == 2:
return constantsimp(Eq(z, C1*zsol[0].rhs + C2*zsol[1].rhs), [C1, C2])
if len(zsol) == 1:
sol1 = zsol[0].rhs
sol2 = sol1*Integral(exp(Integral(-b, x).doit())/sol1**2, x)
zsol.append(Eq(z, sol2))
return constantsimp(Eq(z, C1*zsol[0].rhs + C2*zsol[1].rhs), [C1, C2])
Copy link
Contributor

Choose a reason for hiding this comment

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

We shouldn't use constantsimp here. Also it would be better not to evaluated the integrals automatically or at least it should be possible to control that.

return zsol

def test_kovacic_sol():
f = Function('f')
x = symbols('x')

eq = f(x).diff(x, 2) - 2*f(x)/x**2
sol = find_kovacic_sol(eq)
assert checkodesol(eq, sol)

eq = f(x).diff(x, 2) + (3/(16*x**2*(x - 1)**2))*f(x)
sol = find_kovacic_sol(eq)
assert checkodesol(eq, sol)