From: Tyler Massaro <tyler.massaro_at_gmail.com>

Date: Thu, 24 Jun 2010 13:51:31 -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 Thu 24 Jun 2010 - 18:58:06 GMT

Date: Thu, 24 Jun 2010 13:51:31 -0400

#################

# 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. Received on Thu 24 Jun 2010 - 18:58:06 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 24 Jun 2010 - 23:30:35 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.
*