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: [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.
