From: John Tillinghast <tilling_at_gmail.com>

Date: Wed, 09 Apr 2008 19:18:21 -0400

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. Received on Wed 09 Apr 2008 - 23:21:02 GMT

Date: Wed, 09 Apr 2008 19:18:21 -0400

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. Received on Wed 09 Apr 2008 - 23:21:02 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Thu 10 Apr 2008 - 15:00:28 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*