[R] help, bifurcation diagram efficiency

[Joris Meys]

Basically, don't write loops. Think vectors, matrices,... The R
Inferno of Patrick Burns contains a lot of valuable information on
optimizing code :

http://lib.stat.cmu.edu/S/Spoetry/Tutor/R_inferno.pdf

Cheers
Joris


On Thu, Jun 24, 2010 at 7:51 PM, Tyler Massaro <tyler.massaro at gmail.com> wrote:
> Hello all -
>
> This code will run, but it bogs down my computer when I run it for finer and
> finer time increments and more generations.  I was wondering if there is a
> better way to write my loops so that this wouldn't happen.  Thanks!
>
> -Tyler
>
> #################
>
> # Bifurcation diagram
>
> # Using Braaksma system of equations
>
> # We have however used a Fourier analysis
>
> # to get a forcing function similar to
>
> # cardiac action potential...
>
> #################
>
> require(odesolve)
>
>
>
> # We get this s_of_t function from Maple ws
>
> s_of_t = function(t)
>
> {
>
> (1/10) * (( (1/2) + (1/2) * (sin((1/4)*pi*t))/(abs(sin((1/4)*pi*t)))) * (
> 6.588315815*sin((1/4)*pi*t) - 1.697435362*sin((1/2)*pi*t) - 1.570296922*sin
> ((3/4)*pi*t) + 0.3247901958*sin(pi*t) + 0.7962749105*sin((5/4)*pi*t) +
> 0.07812230515*sin((3/2)*pi*t) - 0.3424877143*sin((7/4)*pi*t) - 0.1148306748*
> sin(2*pi*t) + 0.1063696962*sin((9/4)*pi*t) + 0.02812403009*sin((5/2)*pi*t)))
>
>  }
>
>
> ModBraaksma = function(t, n, p)
>
> {
>
>
>  dx.dt = (1/0.01)*(n[2]-((1/2)*n[1]^2+(1/3)*n[1]^3))
>
>  dy.dt = -(n[1]+p["alpha"]) + 0.032 * s_of_t(t)
>
> list(c(dx.dt, dy.dt))
>
> }
>
>
> initial.values = c(0.1, -0.02)
>
>
> alphamin = 0.01
>
> alphamax = 0.02
>
>
> alphas = seq(alphamin, alphamax, by = 0.00001)
>
>
> TimeInterval = 100
>
>
> times = seq(0.001, TimeInterval, by = 0.001)
>
>
> plot(1, xlim = c(alphamin, alphamax), ylim = c(0,0.3), type = "n",xlab
> = "Values
> of alpha", ylab = "Approximate loop size for a limit cycle", main =
> "Bifurcation
> Diagram")
>
>
> for (i in 1:length(alphas)){
>
>  params = c(alpha=alphas[i])
>
> out = lsoda(initial.values, times, ModBraaksma, params)
>
> X = out[,2]
>
> Y = out[,3]
>
>  for(j in 200:length(times)){
>
>  if (abs(X[j]) < 0.001) {
>
> points(alphas[i], Y[j], pch = ".")
>
>  }
>
> }
>
> }
>
>        [[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.
>



-- 
Joris Meys
Statistical consultant

Ghent University
Faculty of Bioscience Engineering
Department of Applied mathematics, biometrics and process control

tel : +32 9 264 59 87
Joris.Meys at Ugent.be
-------------------------------
Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php