-
-
Notifications
You must be signed in to change notification settings - Fork 38
Incorporate Steven's comments #2
New issue
Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? # to your account
Conversation
gauss([T,] N) | ||
|
||
Return a pair `(x, w)` of `N` quadrature points `x[i]` and weights `w[i]` to | ||
integrate functions on the interval `(-1, 1)`, i.e. `dot(f(x), w)` |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
dot(w, f.(x))
, so that you don't conjugate f
. Or sum(w .* f.(x))
is probably more natural here.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
`T` is an optional parameter specifying the floating-point type, defaulting to `Float64`.
Arbitrary precision (`BigFloat`) is also supported.
Return a pair `(x, w)` of `N` quadrature points `x[i]` and weights `w[i]` to | ||
integrate functions on the interval `(-1, 1)`, i.e. `dot(f(x), w)` | ||
approximates the integral. Uses the method described in Trefethen & | ||
Bau, Numerical Linear Algebra, to find the `N`-point Gaussian quadrature. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
... in O(N
²) operations
Since the rule is symmetric, this only returns the `n+1` points with `x <= 0`. | ||
The function Also computes the embedded `n`-point Gauss quadrature weights `gw` | ||
(again for `x <= 0`), corresponding to the points `x[2:2:end]`. Returns `(x,w,wg)`. | ||
""" |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Given these points and weights, the estimated integral `I` and error `E` can be computed for an
integrand `f(x)` as follows:
x,w,wg = kronrod(n)
fx⁰ = f(x[end]) # f(0)
x⁻ = x[1:end-1] # the x < 0 Kronrod points
fx = f.(x⁻) .+ f.((-).(x⁻)) # f(x < 0) + f(x > 0)
I = sum(fx .* w[1:end-1]) + fx⁰ * w[end]
if isodd(n)
E = abs(sum(fx[2:2:end] .* wg[1:end-1]) + fx⁰*wg[end] - I)
else
E = abs(sum(fx[2:2:end] .* wg[1:end])- I)
end
Also, should have a similar comment to above explaining T
and noting that it is O(N^2).
See also my suggestion on the |
How are things looking now? |
Shouldn't |
Good call, thanks! |
LGTM |
Fantastic. Thanks for all your help! |
I'll put all of the requested changes into this PR so that it can be tagged once this is merged.
Fixes #1