# R/LorenzSystem.R In brasmus/GHE: GHE

```#require(deSolve)

LorenzSystem <- function(X.0=c(-0.1,0.1,0),N=100000,sigma=10,beta=8/3, rho=28) {
# Lorenz force: F = q ( E + v x B)
# Magnetic field flux: M*i_L
# F_Lorentz = M I_L I_B
# F_friction = F omega

dX.dt <- function(t,y,parms) {
with(as.list(parms), {
X <- y[1]; Y <- y[2]; Z <-  y[3]
dx.dt <- sigma*(Y - X)
dy.dt <- X*(rho - Z) - Y
dz.dt <- X*Y - beta*Z
yy <- c(dx.dt,dy.dt,dz.dt)
list(yy)
})
}

y <- X.0
attr(y,'names') <- c("x","y","z")
parms <- c(sigma=sigma, beta=beta, rho=rho)
time <- seq(0,0.001*N,length=N)
out <- as.data.frame(rk4(y,time,dX.dt,parms))

main <- "The Lorenz system"
par(xaxt="n",yaxt="n",bty="n", xpd=NA)
m <- 200
col <- rgb(seq(1,0,length=m)^2,
sin(pi*(1:m)/m)^2,
seq(0,1,length=m)^2)
colind <- round((out\$z - min(out\$z))/(max(out\$z)-min(out\$z))*m)
colind[colind<1] <- 1; colind[colind>m] <- m
colours <- col[colind]

plot(out\$x,out\$y,pch=19,cex=0.3,main=main,
xlab=expression(x),ylab=expression(y),col=colours)

dx <- ( max(out\$x) - min(out\$x) )/10
dy <- ( max(out\$y) - min(out\$y) )/10

legend(max(out\$x)-3*dx,min(out\$y)+dy,
c(expression(sigma==phantom(0)),expression(beta==phantom(0)),
expression(rho==phantom(0)),sigma,beta,rho),bty="n",ncol=2)

text(min(out\$x),max(out\$y),
expression(frac(d*x,dt)==sigma*(y-x)),pos=4)
text(min(out\$x),max(out\$y)-dy,
expression(frac(d*y,dt)==x*(rho -z) -y),pos=4)
text(min(out\$x),max(out\$y)-2*dy,
expression(frac(d*z,dt)==x*y - beta*z),pos=4)
invisible(out)
}

LorenzSystem() -> strangeattractor
```
brasmus/GHE documentation built on May 13, 2019, 2:30 a.m.