[R] LSODA not accurate when RK4 is; what's going on?
Ravi Varadhan
rvaradhan at jhmi.edu
Thu Apr 10 16:13:22 CEST 2008
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.
More information about the R-help
mailing list