John Tillinghast
2008-Apr-09 23:18 UTC
[R] LSODA not accurate when RK4 is; what's going on?
I'm solving the differential equation dy/dx = xy-1 with y(0) = sqrt(pi/2).
This can be used in computing the tail of the normal distribution.
(The actual solution is y(x) = exp(x^2/2) * Integral_x_inf {exp(-t^2/2) dt}
= Integral_0_inf {exp (-xt - t^2/2) dt}. For large x, y ~ 1/x, starting
around x~2.)
I'm testing both lsoda and rk4 from the package odesolve.
rk4 is accurate using step length 10^-2 and probably would be with even
larger steps.
lsoda is pretty accurate out to about x=4, then starts acting strangely. For
step length 10^-3, y suddenly starts to increase after 4, when it should be
strictly decreasing. For step length 10^-4, y instead turns down and start
dropping precipitously.
Any ideas why lsoda would go off the rails when rk4 does so well? I will
soon be using R to solve more complicated systems of ODEs which I don't
understand as well, so I want to know when it can mislead me.
Code:
t4 <- seq(0, 5, by=.0001)> fn
function(t,y,parms=0){return(list(t*y-1))}
s4 <- lsoda(sqrt(pi/2), t4, fn, 0)
[[alternative HTML version deleted]]
John,
You should decrease atol and/or rtol to get accurate integration of your
function.
Try this:
fn <- function(t,y,parms=0){return(list(t*y-1))}
t4 <- seq(0, 5, by=.0001)
s4 <- lsoda(y=sqrt(pi/2), times=t4, func=fn, parms=0, atol=1.e-10,
rtol=1.e-10)
plot(s4, type="l")
Hope this is helpful,
Ravi.
----------------------------------------------------------------------------
-------
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: rvaradhan at jhmi.edu
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
----------------------------------------------------------------------------
--------
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On
Behalf Of John Tillinghast
Sent: Wednesday, April 09, 2008 7:18 PM
To: r-help at r-project.org
Subject: [R] LSODA not accurate when RK4 is; what's going on?
I'm solving the differential equation dy/dx = xy-1 with y(0) = sqrt(pi/2).
This can be used in computing the tail of the normal distribution.
(The actual solution is y(x) = exp(x^2/2) * Integral_x_inf {exp(-t^2/2) dt}
= Integral_0_inf {exp (-xt - t^2/2) dt}. For large x, y ~ 1/x, starting
around x~2.)
I'm testing both lsoda and rk4 from the package odesolve.
rk4 is accurate using step length 10^-2 and probably would be with even
larger steps.
lsoda is pretty accurate out to about x=4, then starts acting strangely. For
step length 10^-3, y suddenly starts to increase after 4, when it should be
strictly decreasing. For step length 10^-4, y instead turns down and start
dropping precipitously.
Any ideas why lsoda would go off the rails when rk4 does so well? I will
soon be using R to solve more complicated systems of ODEs which I don't
understand as well, so I want to know when it can mislead me.
Code:
t4 <- seq(0, 5, by=.0001)> fn
function(t,y,parms=0){return(list(t*y-1))}
s4 <- lsoda(sqrt(pi/2), t4, fn, 0)
[[alternative HTML version deleted]]
______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
John,
I should also mention that by looking at the R source for lsoda() you can
see that the default values for "atol" and "rtol" are
1.e-06. This should
have been mentioned in the help file, but is not.
Another approach to solving your problem is to restrict "hmax", the
maximum
steplength.
s4 <- lsoda(y=sqrt(pi/2), times=t4, func=fn, parms=0, hmax=0.001)
plot(s4, type="l")
Best,
Ravi.
----------------------------------------------------------------------------
-------
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: rvaradhan at jhmi.edu
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
----------------------------------------------------------------------------
--------
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On
Behalf Of John Tillinghast
Sent: Wednesday, April 09, 2008 7:18 PM
To: r-help at r-project.org
Subject: [R] LSODA not accurate when RK4 is; what's going on?
I'm solving the differential equation dy/dx = xy-1 with y(0) = sqrt(pi/2).
This can be used in computing the tail of the normal distribution.
(The actual solution is y(x) = exp(x^2/2) * Integral_x_inf {exp(-t^2/2) dt}
= Integral_0_inf {exp (-xt - t^2/2) dt}. For large x, y ~ 1/x, starting
around x~2.)
I'm testing both lsoda and rk4 from the package odesolve.
rk4 is accurate using step length 10^-2 and probably would be with even
larger steps.
lsoda is pretty accurate out to about x=4, then starts acting strangely. For
step length 10^-3, y suddenly starts to increase after 4, when it should be
strictly decreasing. For step length 10^-4, y instead turns down and start
dropping precipitously.
Any ideas why lsoda would go off the rails when rk4 does so well? I will
soon be using R to solve more complicated systems of ODEs which I don't
understand as well, so I want to know when it can mislead me.
Code:
t4 <- seq(0, 5, by=.0001)> fn
function(t,y,parms=0){return(list(t*y-1))}
s4 <- lsoda(sqrt(pi/2), t4, fn, 0)
[[alternative HTML version deleted]]
______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
John,
In "lsoda" the increment for "time-marching" is adaptively
chosen,
controlled by "atol" and "rtol". The increment is
restricted to lie in
(hmin, hmax). So you can ensure accuracy, perhaps, at the cost of
efficiency by specifying smaller than default values for atol, rtol, and
hmax, which are 1e-06, 1e-06, and Inf, respectively.
"rk4", on the other hand is not so intelligent. It basically use
fixed time
increment.
Ravi.
----------------------------------------------------------------------------
-------
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: rvaradhan at jhmi.edu
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
----------------------------------------------------------------------------
--------
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On
Behalf Of John Tillinghast
Sent: Wednesday, April 09, 2008 7:18 PM
To: r-help at r-project.org
Subject: [R] LSODA not accurate when RK4 is; what's going on?
I'm solving the differential equation dy/dx = xy-1 with y(0) = sqrt(pi/2).
This can be used in computing the tail of the normal distribution.
(The actual solution is y(x) = exp(x^2/2) * Integral_x_inf {exp(-t^2/2) dt}
= Integral_0_inf {exp (-xt - t^2/2) dt}. For large x, y ~ 1/x, starting
around x~2.)
I'm testing both lsoda and rk4 from the package odesolve.
rk4 is accurate using step length 10^-2 and probably would be with even
larger steps.
lsoda is pretty accurate out to about x=4, then starts acting strangely. For
step length 10^-3, y suddenly starts to increase after 4, when it should be
strictly decreasing. For step length 10^-4, y instead turns down and start
dropping precipitously.
Any ideas why lsoda would go off the rails when rk4 does so well? I will
soon be using R to solve more complicated systems of ODEs which I don't
understand as well, so I want to know when it can mislead me.
Code:
t4 <- seq(0, 5, by=.0001)> fn
function(t,y,parms=0){return(list(t*y-1))}
s4 <- lsoda(sqrt(pi/2), t4, fn, 0)
[[alternative HTML version deleted]]
______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.