mathclasstutor

How to solve Lorenz

Attaching package: ‘deSolve’

> Lorenz <- function(t, state, parameters) {
+     with(as.list(c(state, parameters)), {
+         dX <-  a * X + Y * Z
+         dY <-  b * (Y - Z)
+         dZ <- -X * Y + c * Y - Z
+         list(c(dX, dY, dZ))
+     })
+ }
> parameters <- c(a = -8/3, b = -10, c = 28)
> state      <- c(X = 1, Y = 1, Z = 1)
> times      <- seq(0, 100, by = 0.01)
> out <- ode(y = state, times = times, func = Lorenz, parms = parameters)
>
> plot(out)
package ‘scatterplot3d’ successfully unpacked and MD5 sums checked
> library("scatterplot3d", lib.loc="~/R/win-library/3.6")
> scatterplot3d(out[,-1], type = "l")
> logist <- function(t, x, parms) {
+     with(as.list(parms), {
+         dx <- r * x[1] * (1 - x[1]/K)
+         list(dx)
+     })
+ }
> time  <- 0:100
> N0    <- 0.1; r <- 0.5; K <- 100
> parms <- c(r = r, K = K)
> x <- c(N = N0)
> plot(time, K/(1 + (K/N0-1) * exp(-r*time)), ylim = c(0, 120),
+      type = "l", col = "red", lwd = 2)
> time <- seq(0, 100, 2)
> out <- as.data.frame(rk4(x, time, logist, parms))
> points(out$time, out$N, pch = 16, col = "blue", cex = 0.5)
>
> time <- seq(0, 100, 2)
> out <- as.data.frame(euler(x, time, logist, parms))
> points(out$time, out$N, pch = 1)
> time <- seq(0, 100, 4)
> out <- as.data.frame(euler(x, time, logist, parms))
> points(out$time, out$N, pch = 8, cex = 0.5)
>
> out <- as.data.frame(lsoda(x, time, logist, parms))
> points(out$time, out$N, pch = 1, col = "green")
> legend("bottomright",
+        c("analytical","rk4, h=2", "euler, h=2",
+          "euler, h=4", "lsoda"),
+        lty = c(1, NA, NA, NA, NA), lwd = c(2, 1, 1, 1, 1),
+        pch = c(NA, 16, 1, 8, 1),
+        col = c("red", "blue", "black", "black", "green"))
>

> install.packages("marelac", dependencies = FALSE)

package ‘marelac’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
C:\Users\ADMIN\AppData\Local\Temp\Rtmp63NJ0R\downloaded_packages
> air_density(t = 25)
> plot(0:30, air_density(t = 0:30), xlab = "Temperature, dgC", ylab = "kg/m3")
> plot(x= seq(0.8,1.2, 0.01), y = air_density(P = seq(0.8,1.2, 0.01)),
+      xlab = "pressure, bar", ylab = "kg/m3")
 air_density(P = seq(0.8, 1.2, 0.01)) :
  could not find function "air_density"
> image(Bathymetry$x, Bathymetry$y, Bathymetry$z, col = femmecol(100),
+       asp = TRUE, xlab = "", ylab = "") in image(Bathymetry$x, Bathymetry$y, Bathymetry$z, col = femmecol(100),  :
  object 'Bathymetry' not found
>
> air_spechum(t = 25, rh = 50)*1000
Error in air_spechum(t = 25, rh = 50) :
  could not find function "air_spechum"
> atmComp("O2")
Error in atmComp("O2") : could not find function "atmComp"
> atmComp(c("O2", "CH4", "CO2", "N2O"))
Error in atmComp(c("O2", "CH4", "CO2", "N2O")) :
  could not find function "atmComp"
> atmComp()
Error in atmComp() : could not find function "atmComp"
> sum(atmComp())
Error in atmComp() : could not find function "atmComp"
> plot(0:30, viscosity(t = 0:30, S = 35, P = 1),
+      xlab = "temperature", ylab = "g/m/s",
+      main = "shear viscosity of water", type = "l")
> legend("topright", col = c("black","red","blue"), lty = 1,
+        legend = c("S=35, P=1", "S=0, P=1", "S=35, P=100"))
> curve(expr = sin(3*pi*x))
>
> curve(expr = sin(3*pi*x), from = 0, to = 2, col = "blue",
+       xlab = "x", ylab = "f(x)", main = "curve")
> curve(expr = cos(3*pi*x), add = TRUE, col = "red", lty = 2)
>
> abline(h = 0, lty = 2)
>
> legend("bottomleft", c("sin", "cos"),
+        text.col = c("blue", "red"), lty = 1:2)
>
> par(mfrow = c(3, 2))
>
> par(mfcol = c(3, 2))
> par(mfrow = c(2, 2))
>
> for ( i in 1:4) curve(sin(i*pi*x), 0, 1, main = i)