Mercurial > hg > octave-max
changeset 11984:b4986fecdc53 release-3-2-x
Better waypoint transform pour quadgk
| author | Marco Caliari <marco.caliari@univr.it> |
|---|---|
| date | Thu, 11 Jun 2009 07:13:19 +0200 |
| parents | 52b9155fa58a |
| children | bdc383a457fb |
| files | scripts/ChangeLog scripts/general/quadgk.m |
| diffstat | 2 files changed, 9 insertions(+), 7 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,11 +1,15 @@ +2009-06-10 Marco Caliari <marco.caliari@univr.it> + + * general/quadgk.m: Better waypoint transform. + 2009-06-09 David Bateman <dbateman@free.fr> - * general/guadgk.m: Add test case and fixed doubly infinite + * general/quadgk.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 + * general/quadgk.m: Fix doubly infinite transformation to the finite interval. 2009-06-08 Ben Abbott <bpabbott@mac.com>
--- a/scripts/general/quadgk.m +++ b/scripts/general/quadgk.m @@ -193,18 +193,16 @@ ## g(t) = t / (1 - t^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)) + ## t = (2 * g(t)) ./ (1 + sqrt(1 + 4 * g(t) .^ 2)) if (!isempty (waypoints)) - trans = @(x) (-1 + sqrt(1 + 4 * x .^ 2)) ./ (2 * x); + trans = @(x) (2 * x) ./ ((1 + sqrt(1 + 4 * x .^ 2)); 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 (trans(t)) .* (1 + t .^ 2) ./ ((1 - t .^ 2) .^ 2); + f = @(t) f (t ./ (1 - t .^ 2)) .* (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