[R] LSODA not accurate when RK4 is; what's going on?

[Ravi Varadhan]

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.