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

legendre polynomials and gaussian quadrature #313

Merged
merged 22 commits into from
Jun 30, 2021

Conversation

hsnyder
Copy link
Contributor

@hsnyder hsnyder commented Feb 8, 2021

  • Added preliminary implementation for Legendre polynomials
  • Added preliminary implementation for Legendre polynomial first derivatives
  • Added preliminary implementation for Gauss-Legendre quadrature points/weights
  • Added preliminary implementation for Gauss-Legendre-Lobatto quadrature points/weights

link to proposal: #277

- Added preliminary implementation for Legendre polynomial first derivatives
- Added preliminary implementation for Gauss-Legendre quadrature points/weights
- Added preliminary implementation for Gauss-Legendre-Lobatto quadrature points/weights
w(N+1) = 2._dp/(N*(N+1._dp))

do i = 1, int(floor((N+1)/2._dp)-1)
x(i+1) = -cos( (i+0.25_dp)*PI/N - 3/(8*N*PI*(i+0.25_dp)))
Copy link
Member

Choose a reason for hiding this comment

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

A comment containing the source for the initial guess could be helpful.

@certik
Copy link
Member

certik commented Feb 8, 2021

I have a comment regarding the implementation --- what is the accuracy? I had bad experience with higher order (say p=64) that wasn't accurate to machine accuracy with such methods. So we resorted to calculate the weights and points using high precision and generate code that hard codes the numbers valid to all digits. But this comment is an implementation, which is minor.

Let's mainly discuss the API first. I would remove interval as doing [-1, 1] seems like a default?

@ivan-pi
Copy link
Member

ivan-pi commented Feb 8, 2021

Let's mainly discuss the API first. I would remove interval as doing [-1, 1] seems like a default?

The interval is an optional argument:

pure module subroutine gauss_legendre_fp64 (x, w, interval)
        real(dp), intent(out) :: x(:), w(:)
        real(dp), intent(in), optional :: interval(2)

The default is [-1,1]. At least in (chemical) engineering problems, shifted Legendre polynomials are sometimes used. I find the interval argument handy for such cases. Of course one could argue, that the user of the quadrature should be responsible for such things.

I am afraid the intent(out) attribute for x, and w is wrong. If these are allocatable arrays in the calling program, then gauss_legendre will not work as intended. Adding an integer variable, and using fixed size arrays would be the most flexible solution.

@nncarlson
Copy link
Contributor

I am afraid the intent(out) attribute for x, and w is wrong. If these are allocatable arrays in the calling program, then gauss_legendre will not work as intended.

How so? I think you are mistaken here. Whether the actual arguments are allocatable is irrelevant. Are you thinking that the intent(out) causes them to be deallocated?

@ivan-pi
Copy link
Member

ivan-pi commented Feb 8, 2021

@nncarlson, you are right I am mistaken. The concern would only be valid, if the allocatable attribute was given to the dummy variables.

@milancurcic
Copy link
Member

Thank you @hsnyder. I didn't do a proper review yet, but two things stood out on a cursory skim through:

  1. stdlib_functions.f90 for the source file and stdlib_functions for the module are not informative IMO. Is there some descriptive name we can use? What kind of functions are these?
  2. Please use lowercase for all variables. I think uppercase is okay for constants (parameter) although we don't have an established position on this yet.

@ivan-pi
Copy link
Member

ivan-pi commented Feb 8, 2021

  1. Please use lowercase for all variables. I think uppercase is okay for constants (parameter) although we don't have an established position on this yet.

Unfortunately, the lower-case l is sometimes frowned upon due to its potential for mix up with 1. Perhaps p (for the polynomial) and dp (for its derivative) can be used as substitutes.

@awvwgk
Copy link
Member

awvwgk commented Feb 8, 2021

Perhaps p (for the polynomial) and dp (for its derivative) can be used as substitutes.

dp will unfortunately clash with the kind parameter for double precision reals. More than one or two letters might be advisable instead.

@ivan-pi
Copy link
Member

ivan-pi commented Feb 8, 2021

dp will unfortunately clash with the kind parameter for double precision reals. More than one or two letters might be advisable instead.

Right... My next suggestion then are pn and dpn, which are unsurprisingly also used in the code by Stroud & Secrest.

@hsnyder
Copy link
Contributor Author

hsnyder commented Feb 8, 2021

Thanks everyone for the input so far!

Thank you @hsnyder. I didn't do a proper review yet, but two things stood out on a cursory skim through:

1. stdlib_functions.f90 for the source file and `stdlib_functions` for the module are not informative IMO. Is there some descriptive name we can use? What kind of functions are these?

@milancurcic The only thing I can really think of is "special functions". Acceptable?

@arjenmarkus
Copy link
Member

arjenmarkus commented Feb 9, 2021 via email

Copy link
Member

@arjenmarkus arjenmarkus left a comment

Choose a reason for hiding this comment

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

Just a remark: the do-loop at line 33 runs from 0 to newton_iters-1. As the do-loop variable is not explicitly used, I would say: let it run from 1 to newton_iters instead.

@hsnyder
Copy link
Contributor Author

hsnyder commented Feb 10, 2021

Just a remark: the do-loop at line 33 runs from 0 to newton_iters-1. As the do-loop variable is not explicitly used, I would say: let it run from 1 to newton_iters instead.

haha you can tell I'm a C programmer... I will fix this

@hsnyder
Copy link
Contributor Author

hsnyder commented Mar 15, 2021

sorry for the delay I have been very busy. other than resolving the makefile conflicts, what are the next steps towards getting this merged? I will try to get to the conflicts in the makefiles sometime this week but if I can't and someone else wants to, please feel free :)

@hsnyder
Copy link
Contributor Author

hsnyder commented Jun 1, 2021

@awvwgk sorry for the delay. I think test cases and a spec are the outstanding items. In principle anyone could write them but it would obviously be easiest for me to do so. I’ve been very busy for the past few weeks but things are starting to ease up. I think I can finish this off this week. I will post again when I’ve done so.

@hsnyder
Copy link
Contributor Author

hsnyder commented Jun 23, 2021

Sorry for the delay! I've added some tests and specs for the new quadrature functions and the special functions module. Please let me know if any revisions are required.

note: I wasn't sure how to validate the the spec format is correct. I tried to follow the examples already present but I may have made errors and am not sure how to check it.

@awvwgk awvwgk added reviewers needed This patch requires extra eyes and removed waiting for OP This patch requires action from the OP labels Jun 24, 2021
Copy link
Member

@awvwgk awvwgk left a comment

Choose a reason for hiding this comment

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

Thanks a lot, looks good to me now. Only a few comments regarding the markup.

hsnyder and others added 2 commits June 24, 2021 10:20
Copy link
Member

@milancurcic milancurcic left a comment

Choose a reason for hiding this comment

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

Thank you @hsnyder and reviewers!

I have one last minute question and suggestion: What do you think about renaming everything that's called stdlib_specialfunctions* to stdlib_special*? I find specialfunctions quite long and a bit hard to read because two words are together, and special is clear enough what it means IMO.

I think this PR can go forward either way, but I'd prefer if we called this just stdlib_special.

@hsnyder
Copy link
Contributor Author

hsnyder commented Jun 28, 2021

Hi @milancurcic , I prefer specialfunctions because of clarity, but I don't have strong feelings about it. On the one hand, I agree that "special" is shorter, easier to type and read, but on the other hand it's not very descriptive (I don't think "special" on it's own is very clear - the stdlib includes more than just math, so I think "special what?" is a perfectly valid reaction).

I'll wait ~ 24 hours for any additional opinions re: the name change and then I'll merge the request. Thanks very much for reviewing!

@arjenmarkus
Copy link
Member

arjenmarkus commented Jun 29, 2021 via email

@hsnyder
Copy link
Contributor Author

hsnyder commented Jun 29, 2021

So it turns out I cannot merge this myself as I don't have write access. @milancurcic are you able to do so? Feel free to change the name if you want (though I still prefer "specialfunctions", or Arjen's suggestion instead of just "special").

Thanks!

@milancurcic
Copy link
Member

Thanks, @hsnyder, I'll merge as is.

@milancurcic milancurcic merged commit 9fb85ff into fortran-lang:master Jun 30, 2021
@awvwgk awvwgk removed the reviewers needed This patch requires extra eyes label Jul 4, 2021
@mjben
Copy link

mjben commented Jan 13, 2022

Unless I missed something, there should not be x(1) = interval(1) and x(size(x)) = interval(2) in gauss_legendre subroutines, this would only be relevant for gauss_legendre_lobatto.

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.

8 participants