项目作者: shill1729

项目描述 :
numerical solver for first order ODE systems via Euler and RK4 schemes
高级语言: R
项目地址: git://github.com/shill1729/odeSolveR.git
创建时间: 2020-04-24T01:40:47Z
项目社区:https://github.com/shill1729/odeSolveR

开源协议:Other

下载


odeSolveR

A basic ODE solver using Euler scheme and fourth-order Runge-Kutta scheme implementations to numerically solve arbitrary first order one dimensional ODEs and systems of first-order ODEs. The default scheme used is the latter RK4 scheme.

The documentation should be sufficient for setting up custom problems, however, we also provide below a series of examples from the 1970s to the 2000s demonstrating how to use the generic solver. The examples are pulled from S. Strogatz’s excellent Nonlinear Dynamics and Chaos (2015) with the original authors of the model listed.

Table of contents

  1. Installation

  2. One dimensional ODEs

  3. The SIR model

  4. Arbitrary first-order ODE systems

    a. Replicating the hard-coded SIR example

    b. The Lorenz system

    c. Vasquez and Redner political opinion dynamics

    d. Jordan and Smith Keynesian Cross model

    e. Eigen and Schuster Hypercycle equation

    f. Haken two-mode laser model

    g. LRC circuit

  5. SIDARTHE examples

Installation

You can install the latest version of odeSolveR on GitHub with:

  1. devtools::install_github("shill1729/odeSolveR")

One dimensional ODEs

  1. library(odeSolveR)
  2. # The ODE is given in the form y'(t)=f(t, y) with initial y(t_0)=y_0.
  3. # Here we solve y'=cos(t), y(0)=0. The solution is y=sin(t), of course.
  4. f <- function(t, y) cos(t)
  5. y <- ode1(f = f, y0 = 0, tn = 10, n = 5000)
  6. plot(y)

The SIR model

Running an Euler scheme on the SIR model from mathematical epidemiology.

  1. library(odeSolveR)
  2. # Solution to SIR model with 2 days between contacts, 3 days between recovery
  3. # 90 days out,
  4. sir.dat <- sir(beta = 1/2, gamma = 1/3, s0 = 50, i0 = 1, tn = 90, n = 1000)

Arbitrary first-order ODE systems

We can solve arbitrary systems using special syntax. As of now, no error handling is implemented so systems with irregularities, singularities, or blow ups will not be handled with any grace…

Replicating the hard-coded SIR example

This will replicate the above SIR example (without all the plots…) that is implemented as a hard-coded ODE.

  1. library(odeSolveR)
  2. # Numerical parameters
  3. t0 <- 0
  4. tn <- 100
  5. n <- 1000
  6. # Initial conditions
  7. IC <- list(s = 1-0.01, i = 0.01, r = 0)
  8. # Model parameters
  9. parameters <- list(beta = 1/2, gamma = 1/3)
  10. # TODO: need a less hacky fix
  11. list2env(parameters, envir = .GlobalEnv)
  12. # The functions describing rate of change of each state variable in the ODE system
  13. f <- list(function(t, s, i, r) -beta*s*i,
  14. function(t, s, i, r) beta*i*s-gamma*i,
  15. function(t, s, i, r) gamma*i
  16. )
  17. # Solve the system via RK4 iterations
  18. sol <- ode(f, IC, tn = tn, n = n)
  19. # Plotting trajectories
  20. par(mfrow = c(2, 1))
  21. plot_trajectories(sol, legend_names = c("susceptible", "infected", "recovered"))
  22. # Plotting phase portrait
  23. plot_phase_portrait(sol)
  24. plot3D::points3D(gamma/beta, 0, 1, add = TRUE)

SIR

The Lorenz system

As another example of entering in our own ODEs, we can solve the Lorenz system.

  1. library(odeSolveR)
  2. # Numerical parameters
  3. t0 <- 0
  4. tn <- 40
  5. n <- 10000
  6. # Model parameters
  7. parameters <- list(sigma = 10, rho = 28, beta = 8/3)
  8. # TODO: need a less hacky fix to load in parameters so they do not be referenced as parameters$parm in function bodies
  9. list2env(parameters, envir = .GlobalEnv)
  10. # Create the functions that describe the rates of change for each state variable in the ODE system
  11. f <- list(function(t, x, y, z) sigma*(y-x),
  12. function(t, x, y, z) x*(rho-z)-y,
  13. function(t, x, y, z) x*y-beta*z
  14. )
  15. # Initial conditions: names and order must match variables appearing after "t" in function definitions above!
  16. IC <- list(x = 1, y = 1, z = 1)
  17. # Solve the system using RK4
  18. sol <- ode(f, IC, tn = tn, n = n)
  19. # Plot 2 dimension phase space of x-z variables
  20. plot(sol$x, sol$z, type = "l", main = "x-z Phase plane")
  21. # Plot 3 dimensional phase portrait
  22. plot3D::lines3D(x = sol$x, y = sol$y, z = sol$z, main = "3 dimensional phase portrait")

Lorenz

Vasquez and Redner political opinion dynamics

A highly simplified model of political opinion dynamics due to Vasquez and Redner (2004 p. 8489) consisting of a population of rightists (x), leftists (y), and centrists (z). The dynamics folows the assumptions that leftists and rightists never talk to each other but when either talks to a centrist, one of them convinces the other to change his or her political opinion with the winner depending on the sign of a parameter (r). If this parameter is positive, extremists always win and persuade the centrist to move to that end of the spectrum. Otherwise, the centrist pulls the others towards the center.

  1. library(odeSolveR)
  2. # Numerical parameters
  3. t0 <- 0
  4. tn <- 90
  5. n <- 10000
  6. # Model parameters
  7. parameters <- list(r = -0.1)
  8. # TODO: need a less hacky fix to load in parameters so they do not be referenced as parameters$parm in function bodies
  9. list2env(parameters, envir = .GlobalEnv)
  10. # Create the functions that describe the rates of change for each state variable in the ODE system
  11. f <- list(function(t, x, y, z) r*x*z,
  12. function(t, x, y, z) r*y*z,
  13. function(t, x, y, z) -r*x*z-r*y*z
  14. )
  15. # Initial conditions
  16. IC <- list(x = 0.5, y = 0.4, z = 0.10)
  17. # Solve the system using RK4
  18. sol <- ode(f, IC, tn = tn, n = n)
  19. par(mfrow = c(2, 1))
  20. # State trajectories over time
  21. plot(sol$time, sol$x, type = "l", col = "red", ylim = c(0, 1), main = "% of parties over time")
  22. lines(sol$time, sol$y, col = "blue")
  23. lines(sol$time, sol$z, col = "black")
  24. # 3 Dimensional phase portrait
  25. plot3D::lines3D(x = sol$x, y = sol$y, z = sol$z, theta = 120, main = "Phase portrait")

Politics

Jordan and Smith Keynesian Cross model

Adapted from Exercise 2.24 in Jordan and Smith (1987). A simple model for a national economy is given by three variables: income (I), rate of consumer spending (C) and rate of government spending (G), with two parameters alpha and beta both greater than unity. The parameters here correspond to a semi-oscillating trajectory that takes a long time to reach equillibrium.

  1. library(odeSolveR)
  2. # Numerical parameters
  3. t0 <- 0
  4. tn <- 300
  5. n <- 5000
  6. # Model parameters
  7. parameters <- list(alpha = 1.01, beta = 1.01, k = 0.2)
  8. # TODO: need a less hacky fix to load in parameters so they do not be referenced as parameters$parm in function bodies
  9. list2env(parameters, envir = .GlobalEnv)
  10. # Create the functions that describe the rates of change for each state variable in the ODE system
  11. f <- list(function(t, I, C, G) I-alpha*C,
  12. function(t, I, C, G) beta*(I-C-G),
  13. function(t, I, C, G) 0 # Or use k*(I-alpha*C) for G=G0+k*I linear growth of expenditures in terms of income
  14. # or use 2*k*I*(I-alpha*C) for quadratic growth of expenditures in terms of income
  15. )
  16. # Initial conditions
  17. IC <- list(I = 300, C = 100, G = 300*k)
  18. # Solve the system using RK4
  19. sol <- ode(f, IC, tn = tn, n = n)
  20. par(mfrow = c(2, 1))
  21. # State trajectories over time
  22. plot(sol$time, sol$I, type = "l", col = "blue")
  23. lines(sol$time, sol$C, col = "red")
  24. lines(sol$time, sol$G, col = "green")
  25. # Plot 2 dimension phase space of x-z variables
  26. plot(sol$I, sol$C, type = "l")

Economy

Eigen and Schuster Hypercycle equation

In Eigen and Schuster’s (1978) model of pre-biotic evolution, a group of RNA molecules or other self-replicating chemical units are imagined to catalyze each other’s reproduction in a closed feedback loop, with one molecule serving as the catalyst for the next. The simplest scheme considered by Eigen and Schuster’s is implemented here, in dimensionless form as the hypercycle equations, for three molecules. For the dynamics of the hypercycle equation under larger values of the number of molecules we refer interested readers to Hofbauer and Sigmund (1998).

  1. library(odeSolveR)
  2. # Numerical parameters
  3. t0 <- 0
  4. tn <- 45
  5. n <- 5000
  6. parameters <- list()
  7. # Create the functions that describe the rates of change for each state variable in the ODE system
  8. f <- list(function(t, x1, x2, x3) x1*(x3-x1*x3-x2*x1-x3*x2),
  9. function(t, x1, x2, x3) x2*(x1-x1*x3-x2*x1-x3*x2),
  10. function(t, x1, x2, x3) x3*(x2-x1*x3-x2*x1-x3*x2)
  11. )
  12. # Initial conditions
  13. IC <- list(x1 = 0.3, x2 = 0.6, x3 = 0.1)
  14. # Solve the system using RK4
  15. sol <- ode(f, IC, tn = tn, n = n)
  16. par(mfrow = c(2, 1))
  17. # State trajectories over time
  18. plot(sol$time, sol$x1, type = "l", col = "blue", ylim = c(0, 1), main = "Relative frequency of molecules")
  19. lines(sol$time, sol$x2, col = "red")
  20. lines(sol$time, sol$x3, col = "black")
  21. # 3 Dimensional phase portrait
  22. plot3D::lines3D(x = sol$x1, y = sol$x2, z = sol$x3, main = "Phase portrait")

Hypercycle

Haken two-mode laser model

According to Haken (1983, p. 129), a two-mode laser produces two different kinds of photons. The number of each kind of photon is modeled via a first order ODE system, implemented below. The rate of change of the number of photons takes a basic “gain-loss” form dependent on the number of excited photons.

  1. library(odeSolveR)
  2. # Numerical parameters
  3. t0 <- 0
  4. tn <- 45
  5. n <- 5000
  6. # Model parameters
  7. parameters <- list(G1 = 0.5, G2 = 0.4, k1 = 0.3, k2 = 0.2, alpha1 = 0.9, alpha2 = 0.02, N0 = 100)
  8. # TODO: need a less hacky fix to load in parameters so they do not be referenced as parameters$parm in function bodies
  9. list2env(parameters, envir = .GlobalEnv)
  10. # Create the functions that describe the rates of change for each state variable in the ODE system
  11. f <- list(function(t, n1, n2) G1*(N0-alpha1*n1-alpha2*n2)-k1*n1,
  12. function(t, n1, n2) G2*(N0-alpha1*n1-alpha2*n2)-k2*n2
  13. )
  14. # Initial conditions
  15. IC <- list(n1 = 0, n2 = 0)
  16. # Solve the system using RK4
  17. sol <- ode(f, IC, tn = tn, n = n)
  18. par(mfrow = c(2, 1))
  19. # State trajectories over time
  20. plot(sol$time, sol$n1, type = "l", col = "blue", ylim = c(min(sol$n1, sol$n2), max(sol$n1, sol$n2)), main = "Number of photons")
  21. lines(sol$time, sol$n2, col = "red")
  22. # 2 Dimensional phase portrait
  23. plot(sol$n1, sol$n2, main = "Phase portrait", type = "l")

Laser

LRC circuit

The charge of the capacitor in a basic LRC circuit can be modeled by a second-order ODE. Here the (x) coordinate is the charge, and (y) is the current (the time derivative of the charge.)

  1. library(odeSolveR)
  2. tn <- 20
  3. n <- 4000
  4. # Initial conditions
  5. #' x = charge of the capacitor
  6. #' y = current (time-derivative of charge, i.e. velocity of charge)
  7. IC <- list(x = 0, y = 0)
  8. # Parameters
  9. #' L = inductance in henrys
  10. #' R = resistance in ohms
  11. #' C = capicitance in farads
  12. parameters <- list(L = 1,
  13. R = 1,
  14. C = 1/30
  15. )
  16. # Load parameters to global environment
  17. list2env(x = parameters, envir = .GlobalEnv)
  18. # Energy, the non-homogeneous component of the ODE
  19. energy <- function(t) {100}
  20. # Standard form functions of ODE
  21. f <- list(
  22. function(t, x, y) y,
  23. function(t, x, y) energy(t)-(R/L)*y-x/(C*L)
  24. )
  25. # RK4 Solution
  26. sol <- ode(f, IC = IC, t0 = 0, tn = tn, n = n)
  27. # Summary of specs:
  28. # in units of angular freq. Measures how fast the transient response of the circuit will die away after stimulus is removed
  29. attenuation <- R/(2*L) # nepers per second
  30. angular_resonance_freq <- 1/sqrt(L*C) # angular frequency
  31. damping_factor <- attenuation/angular_resonance_freq # unitless damping factor
  32. specs <- data.frame(inductance = L, # henry
  33. resistance = R, # Ohms
  34. capacitance = C, # farads
  35. attenuation = attenuation, # Nepers per second
  36. bandwidth = 2*attenuation, # nepers per second
  37. fractional_bandwith = 2*attenuation/angular_resonance_freq, # percentage
  38. angular_resonance_freq = angular_resonance_freq, # angular frequency
  39. resonance_freq = angular_resonance_freq/(2*pi), # hertz
  40. damping_factor = damping_factor # unitless
  41. )
  42. specs$transient_type <- ifelse(damping_factor > 1, "overdamped", ifelse(damping_factor < 1, "underdamped", "critically damped"))
  43. specs$Q_factor <- 1/specs$fractional_bandwith
  44. specs$damped_resonance_freq <- ifelse(damping_factor < 1, sqrt(angular_resonance_freq^2-attenuation^2), NA)
  45. # Linear- matrix and stability classification
  46. fixed_point <- c(C*L*energy(sol$time[n+1]), 0)
  47. A <- rbind(c(0, 1), c(-1/(C*L), -R/L))
  48. # Output:
  49. # Plotting trajectories and phase-plane
  50. par(mfrow = c(2, 1))
  51. plot_trajectories(sol, legend_loc = "topright", legend_names = c("charge", "current"), legend_size = 0.7)
  52. plot_phase_portrait(sol)
  53. points(sol$x[1], sol$y[1], col = "green")
  54. points(C*L*energy(sol$time[n+1]), 0, col = "red")
  55. # Final print
  56. print(specs)
  57. print(classify_equilibrium(A))
  58. print(fixed_point)

LRC circuit