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.mass...@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@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.m...@ugent.be ------------------------------- Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php ______________________________________________ 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.
