Skip to content

More Gaussian Quadrature Methods - #1517

Open
JacobHass8 wants to merge 37 commits into
boostorg:developfrom
JacobHass8:orthogonal-polynomial
Open

JacobHass8 wants to merge 37 commits into
boostorg:developfrom
JacobHass8:orthogonal-polynomial

Conversation

@JacobHass8

Copy link
Copy Markdown
Contributor

I'm working on adding more Gaussian quadrature methods using Laguerre, Associated Laguerre and Jacobi polynomials. This is done by generalizing the orthogonal polynomial root finder in #1461 to any family of polynomials. I also added tests for finding the root of these polynomials in the process.

So far I've only implemented everything for Laguerre polynomials, but I'm going to start on Associate Laguerre polynomials next.

@NAThompson

Copy link
Copy Markdown
Collaborator

@JacobHass8 : I encourage you to take on the adaptive versions of these in the first pass-or at least give an error estimate. This truly increases the value of these methods to users.

@JacobHass8

Copy link
Copy Markdown
Contributor Author

@JacobHass8 : I encourage you to take on the adaptive versions of these in the first pass-or at least give an error estimate. This truly increases the value of these methods to users.

I deleted my earlier comment because I realized adaptive methods are much, much harder for Hermite and Laguerre polynomials than simply using Gauss-Kronrod. I went back to your previous comment and took a look at the paper here: https://arxiv.org/html/1203.1554v2. Could I just use the Mathematica code provided to generate the abscissa and weights for fixed number of points? I can also work with Claude to translate the code but I'm not sure how long that would take. It also needs to solve a linear system of equations so will require Eigen or a linear solver.

As you noted earlier, an error estimate could just be provided by evaluating a higher order method and taking the difference. This is very minimal but might not be the best estimation.

Do you have any preference on which route to go?

@NAThompson

Copy link
Copy Markdown
Collaborator

@JacobHass8 : Hmm, I think you're right-the literature is considerably less developed for this problem.

What about tackling Genz-Kaister first?

@JacobHass8

Copy link
Copy Markdown
Contributor Author

@JacobHass8 : Hmm, I think you're right-the literature is considerably less developed for this problem.

What about tackling Genz-Kaister first?

It looks like Genz-Keister will be difficult/impossible to compute on the fly. Apparently the linear system of equations can be ill-conditioned and return complex roots. I think that I'll just use the mathematica code in the arxiv paper to calculate the coefficients to higher precision

@NAThompson

Copy link
Copy Markdown
Collaborator

@JacobHass8 : Could you get your quadrature methods added to the nanobind layer so I can play around with them?

@JacobHass8

Copy link
Copy Markdown
Contributor Author

@JacobHass8 : Could you get your quadrature methods added to the nanobind layer so I can play around with them?

Yes, I can work on that next. I've only really implemented Gauss-Laguerre though. I still need to do Associated Laguerre and Jacobi methods.

I also have weights and abscissa for the nested hermite method. I still need to test if they actually give higher accuracy than regular gauss hermite. I assumed this would go in a separate PR but I can commit it here if you want.

@NAThompson

Copy link
Copy Markdown
Collaborator

I assumed this would go in a separate PR but I can commit it here if you want.

Separate PR is fine.

Comment thread test/test_gauss_laguerre_quadrature_eigen.cpp Outdated
@JacobHass8

Copy link
Copy Markdown
Contributor Author

@NAThompson this is ready for review now. I need to update the documentation but I can do that quickly tonight. I also need to condense some of the code. I copy/pasted the gauss_hermite class to use for gauss_jacobi and gauss_laguerre so there's duplicate code.

One implementation choice I'd like to point out. Gauss-Jacobi requires N, alpha and beta to define the integration rule. Since alpha and beta can be any number greater than -1, I didn't save any nodes/weights. They are all generated on the fly. Thus, the integrate method in the gauss_jacobi class isn't static. This breaks the Gauss quadrature design, but I think that it's warranted here.

@JacobHass8

Copy link
Copy Markdown
Contributor Author

Ah, I still have to add the nanobind layer!

@NAThompson

Copy link
Copy Markdown
Collaborator

@JacobHass8 :

Ah, I still have to add the nanobind layer!

I can just have claude generate one-would you mind that?

@JacobHass8

Copy link
Copy Markdown
Contributor Author

I can just have claude generate one-would you mind that?

That would be greatly appreciated. It would save the headache of figuring out nanobind although it looks similar to pybind11

@NAThompson

Copy link
Copy Markdown
Collaborator

@JacobHass8 : This is next on my plate. I'll begin to mount an attack on it soon.

@NAThompson

Copy link
Copy Markdown
Collaborator

@JacobHass8 : I would like to extend this to the anti-Gaussian quadrature formulas-would you mind me making that extension to your PR?

@JacobHass8

Copy link
Copy Markdown
Contributor Author

@JacobHass8 : I would like to extend this to the anti-Gaussian quadrature formulas-would you mind me making that extension to your PR?

Wow, that's pretty cool! Go for it. Hopefully the orthogonal_polynomial class I made can be adapted to this.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants