Refactor bessel solver by Mohitbalwani26 · Pull Request #21706 · sympy/sympy

Mohitbalwani26 · 2021-07-05T10:12:08Z

References to other Issues or PRs

Brief description of what is fixed or changed

SEE #18348

Other comments

Release Notes

NO ENTRY

sympy-bot · 2021-07-05T10:12:09Z

✅

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.

No release notes entry will be added for this pull request.

Click here to see the pull request description that was parsed.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->


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

#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
NO ENTRY
<!-- END RELEASE NOTES -->

Mohitbalwani26 · 2021-07-05T15:34:55Z

+                    self.rn = {'b':self.rn[a4],'m':self.rn[b4]}
+                    does_match = True
+        return does_match
+


@oscarbenjamin According to my understanding about 2nd_linear_airy, I attempted to simplify the match function something like :

def _matches(self): eq = self.ode_problem.eq_high_order_free f = self.ode_problem.func order = self.ode_problem.order x = self.ode_problem.sym df = f.diff(x) a4 = Wild('a4', exclude=[x,f,df]) b4 = Wild('b4', exclude=[x,f,df]) match = self.ode_problem.get_linear_coefficients(eq, f, order) does_match = False if order == 2 and match and match[2] != 0: if match[1].is_zero: self.rn = match[0].match(a4+b4*x) if self.rn and self.rn[b4] != 0: self.rn = {'b':self.rn[a4],'m':self.rn[b4]} does_match = True return does_match

It passes all the tests but wanted to ask you does this diff look ok?

Does get_linear_coefficients divide by the leading coefficient?

What happens if you have something like x y'' + x y = 0?

x y'' + x y = 0

It fails for this. I thought self.ode_problem.eq_high_order_free passing in the get_linear_coefficients would work as eq_high_order_free removes leading coefficient but it only remove f(x) not x.

if we replace self.rn = match[0].match(a4+b4*x) with self.rn = cancel(match[0]/match[2]).match(a4+b4*x) . It will solve the problem 🤔

asmeurer · 2021-07-06T17:42:31Z

+    def _get_general_solution(self, *, simplify_flag: bool = True):
+        f = self.ode_problem.func.func
+        x = self.ode_problem.sym
+        (C0, C1) = self.ode_problem.get_numbered_constants(num=2)


Shouldn't this be C1 and C2?

Mohitbalwani26 · 2021-07-06T19:19:16Z

@oscarbenjamin I have updated the match function. PTAL

oscarbenjamin · 2021-07-06T20:05:35Z

I just noticed this bug on master:

In [43]: dsolve(x*f(x).diff(x, 2) - x*f(x))
Out[43]: f(x) = √x⋅(C₁⋅besselj(1/2, ⅈ⋅x) + C₂⋅bessely(1/2, ⅈ⋅x))

Is that bug also in this PR?

Mohitbalwani26 · 2021-07-06T20:43:54Z

dsolve(x*f(x).diff(x, 2) - x*f(x))

Isn't it the correct solution as checkodesol(eq,sol) also returns [(True, 0)]?
According to the reference mentioned in the docstring Link it is reducible to bessel equation (4th point)
x²y'' + axy' + (x² - n²) can be reduced to bessel form by substituting y = x^(1-a)/2 z

oscarbenjamin · 2021-07-06T21:14:57Z

Isn't it the correct solution as checkodesol(eq,sol) also returns [(True, 0)]?

It might be correct for some part of the domain but it is definitely not as good as this:

In [68]: eq = f(x).diff(x, 2) - f(x)

In [69]: dsolve(eq)
Out[69]: 
           -x       x
f(x) = C₁⋅ℯ   + C₂⋅ℯ 

In [70]: dsolve(expand(x*eq))
Out[70]: f(x) = √x⋅(C₁⋅besselj(1/2, ⅈ⋅x) + C₂⋅bessely(1/2, ⅈ⋅x))

The Bessel solution is not continuous at x=0 and is significantly more complicated then the basic exponential solution here.

Mohitbalwani26 · 2021-07-06T21:24:53Z

The Bessel solution is not continuous at x=0 and is significantly more complicated then the basic exponential solution here.

So should we first check if the equation is analytic at 0 or some point? 🤔

oscarbenjamin · 2021-07-06T21:26:06Z

So should we first check if the equation is analytic at 0 or some point? 🤔

This should definitely match the constant coefficients solver. I'm not sure what should be done with the Bessel solver.

Mohitbalwani26 · 2021-07-06T21:39:13Z

This should definitely match the constant coefficients solver. I'm not sure what should be done with the Bessel solver.

It doesn't not match directly with constant coefficients solver because in eq_higher_order_free, we are removing the leading coefficient having f(x). Although this equation has x as common so it should match with Factorable solver and then recursively call the constant coefficient solver.

Mohitbalwani26 · 2021-07-07T08:58:48Z

It doesn't not match directly with constant coefficients solver because in eq_higher_order_free, we are removing the leading coefficient having f(x). Although this equation has x as common so it should match with Factorable solver and then recursively call the constant coefficient solver.

@oscarbenjamin Thanks for pointing out this example!
This was a bug and got introduced while refactoring the Factorable hint. As earlier in sympy version 1.5 it was classified by factorable hint and solution given was correct.
I have fixed it now and added the example.

oscarbenjamin · 2021-07-07T12:28:11Z



+def test_slow_examples_factorable():
+    _ode_solver_test(_get_examples_ode_sol_factorable, run_slow_test=True)


Shouldn't this be marked as slow?

oscarbenjamin · 2021-07-07T18:28:29Z

+    Gives solution of the Bessel differential equation
+
+    .. math :: x^2 \frac{d^2y}{dx^2} + x \frac{dy}{dx} y(x) + (x^2-n^2) y(x)
+
+    if n is integer then the solution is of the form Eq(f(x), C0 besselj(n,x)
+    + C1 bessely(n,x)) as both the solutions are linearly independent else if
+    n is a fraction then the solution is of the form Eq(f(x), C0 besselj(n,x)
+    + C1 besselj(-n,x)) which can also transform into Eq(f(x), C0 besselj(n,x)
+    + C1 bessely(n,x)).


The equations here should be in LaTeX. If they are supposed to be code rather than equations then they should be in double backticks.

oscarbenjamin · 2021-07-07T18:34:03Z

        NthLinearEulerEqNonhomogeneousVariationOfParameters: ('nth_linear_euler_eq_nonhomogeneous_variation_of_parameters',),
        NthLinearEulerEqNonhomogeneousUndeterminedCoefficients: ('nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients',),
+        LinearBessel2nd: ('2nd_linear_bessel',),
+        LinearAiry2nd: ('2nd_linear_airy',),


It would be better to call these e.g. SecondLinearBessel

Also Hypergeometric2nd should be renamed.

oscarbenjamin · 2021-07-08T10:48:35Z

Looks good. Thanks

Mohitbalwani26 · 2021-07-08T10:52:11Z

@oscarbenjamin should I raise PR for lie_group first and then speedup the process of calling classify_ode?

oscarbenjamin · 2021-07-08T11:33:52Z

@oscarbenjamin should I raise PR for lie_group first and then speedup the process of calling classify_ode?

Up to you. You can also just work on both in parallel.

refactor bessel solver

86854cc

Mohitbalwani26 requested a review from oscarbenjamin as a code owner July 5, 2021 10:12

refactor airy solver

b265bc8

Mohitbalwani26 commented Jul 5, 2021

View reviewed changes

asmeurer reviewed Jul 6, 2021

View reviewed changes

update match function of airy

96a2654

Mohitbalwani26 added 2 commits July 7, 2021 14:19

fix bug in factorable hint

7688f4e

add example

85e6339

oscarbenjamin reviewed Jul 7, 2021

View reviewed changes

mark slow test function slow

df9cc4f

Mohitbalwani26 requested a review from oscarbenjamin July 7, 2021 17:28

oscarbenjamin reviewed Jul 7, 2021

View reviewed changes

fix docstring and change class name

845bfcb

oscarbenjamin added the solvers.dsolve label Jul 7, 2021

rename hypergeometric class

62b7910

Mohitbalwani26 requested a review from oscarbenjamin July 8, 2021 00:17

oscarbenjamin merged commit 17f3457 into sympy:master Jul 8, 2021

Mohitbalwani26 deleted the refactor_bessel_airy branch July 13, 2021 15:35



		def test_slow_examples_factorable():
		_ode_solver_test(_get_examples_ode_sol_factorable, run_slow_test=True)

Uh oh!

Uh oh!

Conversation

Mohitbalwani26 commented Jul 5, 2021

References to other Issues or PRs

Brief description of what is fixed or changed

Other comments

Release Notes

Uh oh!

sympy-bot commented Jul 5, 2021

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Mohitbalwani26 Jul 5, 2021 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Mohitbalwani26 commented Jul 6, 2021

Uh oh!

oscarbenjamin commented Jul 6, 2021

Uh oh!

Mohitbalwani26 commented Jul 6, 2021

Uh oh!

oscarbenjamin commented Jul 6, 2021

Uh oh!

Mohitbalwani26 commented Jul 6, 2021

Uh oh!

oscarbenjamin commented Jul 6, 2021

Uh oh!

Mohitbalwani26 commented Jul 6, 2021

Uh oh!

Mohitbalwani26 commented Jul 7, 2021

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

oscarbenjamin commented Jul 8, 2021

Uh oh!

Mohitbalwani26 commented Jul 8, 2021

Uh oh!

oscarbenjamin commented Jul 8, 2021

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

Mohitbalwani26 Jul 5, 2021 •

edited

Loading