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: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html ---------------------------------------------------------------------------- -------- -----Original Message----- From: [EMAIL PROTECTED] [mailto:[EMAIL PROTECTED] On Behalf Of John Tillinghast Sent: Wednesday, April 09, 2008 7:18 PM To: r-help@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@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. ______________________________________________ R-help@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.
