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.

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

[GSoC] Kovacic Solver - Part 1 #21770

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

Conversation

naveensaigit
Copy link
Contributor

This PR adds the part 1 of 3 of the Kovacic algorithm, an algorithm for solving second order linear homogeneous ODEs.

References to other Issues or PRs

Brief description of what is fixed or changed

Other comments

Release Notes

  • solvers
    • Part 1 of 3 for Kovacic's algorithm, an algorithm for solving second order linear homogeneous ODEs has been added.

@sympy-bot
Copy link

Hi, I am the SymPy bot (v161). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

  • solvers
    • Part 1 of 3 for Kovacic's algorithm, an algorithm for solving second order linear homogeneous ODEs has been added. (#21770 by @naveensaigit)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.9.

Click here to see the pull request description that was parsed.
This PR adds the part 1 of 3 of the Kovacic algorithm, an algorithm for solving second order linear homogeneous ODEs.

#### References to other Issues or PRs

#### Brief description of what is fixed or changed

#### Other comments

#### Release Notes
<!-- BEGIN RELEASE NOTES -->
* solvers
  * Part 1 of 3 for Kovacic's algorithm, an algorithm for solving second order linear homogeneous ODEs has been added.
<!-- END RELEASE NOTES -->

@sympy-bot
Copy link

🟠

Hi, I am the SymPy bot (v161). I've noticed that some of your commits add or delete files. Since this is sometimes done unintentionally, I wanted to alert you about it.

This is an experimental feature of SymPy Bot. If you have any feedback on it, please comment at sympy/sympy-bot#75.

The following commits add new files:

  • e12d3f4:
    • sympy/solvers/ode/kovacic.py

If these files were added/deleted on purpose, you can ignore this message.


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?

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?

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.

@oscarbenjamin
Copy link
Contributor

We should add a function for reduction of order:
https://en.wikipedia.org/wiki/Reduction_of_order

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

5 participants