Mercurial > hg > octave-nkf
changeset 9323:8a3afa8ebb06
Fixes for quadgk for doubly infinite interval
author | David Bateman <dbateman@free.fr> |
---|---|
date | Tue, 09 Jun 2009 23:18:30 +0200 |
parents | 66fdc831c580 |
children | 29563379fa9b |
files | scripts/ChangeLog scripts/general/quadgk.m |
diffstat | 2 files changed, 17 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,13 @@ +2009-06-09 David Bateman <dbateman@free.fr> + + * general/guadgk.m: Add test case and fixed doubly infinite + waypoint transform for x = 0 case. + +2009-06-09 Marco Caliari <marco.caliari@univr.it> + + * general/guadgk.m: Fix doubly infinite transformation to the finite + interval. + 2009-06-08 Ben Abbott <bpabbott@mac.com> * plot/axis.m: Fix bug for 'axis tight' with multiple lines, modify
--- a/scripts/general/quadgk.m +++ b/scripts/general/quadgk.m @@ -77,7 +77,7 @@ ## function to integrate, then these can be flagged with the ## @code{"WayPoints"} property. This forces the ends of a sub-interval ## to fall on the breakpoints of the function and can result in -## significantly improved estimated of the error in the integral, faster +## significantly improved estimation of the error in the integral, faster ## computation or both. For example, ## ## @example @@ -191,18 +191,20 @@ ## \int_{-\infinity_^\infinity f(x) dx = \int_-1^1 f (g(t)) g'(t) dt ## where ## g(t) = t / (1 - t^2) - ## g'(t) = 1 / (1 + t^2) ^ 2 + ## g'(t) = (1 + t^2) / (1 - t^2) ^ 2 ## waypoint transform is then ## t = (-1 + sqrt(1 + 4 * g(t) .^ 2)) ./ (2 * g(t)) if (!isempty (waypoints)) trans = @(x) (-1 + sqrt(1 + 4 * x .^ 2)) ./ (2 * x); subs = [-1; trans(waypoints); 1]; + ## The waypoint transform gives NaN for x=0; + subs(isnan(subs)) = 0; else subs = linspace (-1, 1, 11)'; endif h = 2; trans = @(t) t ./ (1 - t.^2); - f = @(t) f (t ./ (1 + t .^ 2)) ./ ((1 + t .^ 2) .^ 2); + f = @(t) f (trans(t)) .* (1 + t .^ 2) ./ ((1 - t .^ 2) .^ 2); elseif (isinf(a)) ## Formula defined in Shampine paper as two separate steps. One to ## weaken singularity at finite end, then a second to transform to @@ -442,3 +444,5 @@ %!assert (quadgk (@(x) abs (1 - x.^2), 0, 2, 'Waypoints', 1), 2, 1e-6) %!assert (quadgk(@(x) 1./(sqrt(x).*(x+1)), 0, Inf), pi, 1e-6) %!assert (quadgk (@(z) log (z), 1+1i, 1+1i, 'WayPoints', [1-1i, -1,-1i, -1+1i]), -pi * 1i, 1e-6) + +%!assert (quadgk (@(x) exp(-x .^ 2), -Inf, Inf), sqrt(pi), 1e-6) \ No newline at end of file