Comprehensive Guide to Numerical Integration and ODEs in R

Mastering deSolve, odeintr, and More for Complex ODE Integration

R Progaming
Mathematics
Author

Abhirup Moitra

Published

June 12, 2024

Solving Stiff Differential Equations

deSolve() Package

Description: The deSolve package is a comprehensive tool for solving differential equations in R. It is capable of handling a variety of differential equations, including ordinary differential equations (ODEs), partial differential equations (PDEs), and delay differential equations (DDEs). The package includes several solvers designed for stiff systems, such as lsoda, lsode, lsodes, vode, and daspk.

Example:

library(deSolve)

# Define the stiff ODE function
stiffODE <- function(t, y, parms) {
  dy1 <- -1000 * y[1] + 3000 - 2000 * y[2]
  dy2 <- y[1] - y[2] - 0.5 * y[2]^3
  list(c(dy1, dy2))
}

# Initial state and parameters
state <- c(y1 = 1, y2 = 0)
times <- seq(0, 1, by = 0.01)

# Solve the stiff ODE using lsoda solver
out <- ode(y = state, times = times, func = stiffODE, parms = NULL, method = "lsoda")
plot(out)

odeintr() Package

Description: The odeintr package is designed for integrating differential equations in R, utilizing the Boost.odeint package as its integration engine. This package offers several notable features, including simple specification of the ODE system and named, dynamically settable system parameters. It provides intelligent defaults, which are easily overridden, and supports a wide range of integration methods for compiled systems through various stepper types.

One of the standout features of odeintr is its fully automated compilation of ODE systems specified in C++. It also supports simple OpenMP vectorization for large systems, although this feature is currently broken in the latest Boost release. Results from integrations are returned as a straightforward data frame, ready for analysis and plotting.

The package allows users to specify custom observers in R, which can return arbitrary data. There are three options for calling the observer: at regular intervals, after each update step, or at specified times. Users can also alter the system state and restart simulations from where they left off. Additionally, odeintr can compile an implicit solver with symbolic evaluation of the Jacobian, and it provides the ability to save and edit the generated C++ code easily.

#installation 
#install.packages(odeintr)                   # released
# Install and load the package
devtools::install_github("thk686/odeintr", force = TRUE)
Using GitHub PAT from the git credential store.
Downloading GitHub repo thk686/odeintr@HEAD

── R CMD build ─────────────────────────────────────────────────────────────────
* checking for file 'C:\Users\maths\AppData\Local\Temp\Rtmpm41qZK\remotes1a2046713a2a\thk686-odeintr-1f1e7f9/DESCRIPTION' ... OK
* preparing 'odeintr':
* checking DESCRIPTION meta-information ... OK
* cleaning src
* checking for LF line-endings in source and make files and shell scripts
* checking for empty or unneeded directories
* building 'odeintr_1.7.1.tar.gz'
Warning: file 'odeintr/configure' did not have execute permissions: corrected

Installing package into 'C:/Users/maths/AppData/Local/R/win-library/4.3'
(as 'lib' is unspecified)
library(odeintr)

# Define the ODE function
dxdt = function(x, t) x * (1 - x)

# Integrate the system
system.time({
  x = integrate_sys(dxdt, 0.001, 15, 0.01)
})
   user  system elapsed 
   0.02    0.00    0.02 
# Plot the results
plot(x, type = "l", lwd = 3, col = "steelblue", main = "Logistic Growth")

# Compile the system
compile_sys("logistic", "dxdt = x * (1 - x)")

# Integrate the compiled system
system.time({
  x = logistic(0.001, 15, 0.01)
})
   user  system elapsed 
      0       0       0 
# Plot the results of the compiled system
plot(x, type = "l", lwd = 3, col = "steelblue", main = "Logistic Growth")
points(logistic_at(0.001, sort(runif(10, 0, 15)), 0.01), col = "darkred")

dxdt = function(x, t) c(x[1] - x[1] * x[2], x[1] * x[2] - x[2])
obs = function(x, t) c(Prey = x[1], Predator = x[2], Ratio = x[1] / x[2])
system.time({x = integrate_sys(dxdt, rep(2, 2), 20, 0.01, observer = obs)})
   user  system elapsed 
   0.03    0.00    0.03 
plot(x[, c(2, 3)], type = "l", lwd = 2, col = "steelblue", main = "Lotka-Volterra Phase Plot")

plot(x[, c(1, 4)], type = "l", lwd = 2, col = "steelblue", main = "Prey-Predator Ratio")

# C++ code
Lorenz.sys = '
  dxdt[0] = 10.0 * (x[1] - x[0]);
  dxdt[1] = 28.0 * x[0] - x[1] - x[0] * x[2];
  dxdt[2] = -8.0 / 3.0 * x[2] + x[0] * x[1];
  ' # Lorenz.sys
compile_sys("lorenz", Lorenz.sys)
system.time({x = lorenz(rep(1, 3), 100, 0.001)})
   user  system elapsed 
   0.04    0.00    0.03 
plot(x[, c(2, 4)], type = 'l', col = "steelblue", main = "Lorenz Attractor")

VdP.sys = '
dxdt[0] = x[1];
dxdt[1] = 2.0 * (1 - x[0] * x[0]) * x[1] - x[0];
' # Vdp.sys
compile_sys("vanderpol", VdP.sys, method = "bsd") # Bulirsch-Stoer
system.time({x = vanderpol(rep(1e-4, 2), 100, 0.01)})
   user  system elapsed 
      0       0       0 
par(mfrow = c(2, 2), mar = rep(0.5, 4), oma = rep(5, 4), xpd = NA)
make.plot = function(xy, xlab = NA, ylab = NA)
  plot(xy, col = "steelblue", lwd = 2, type = "l",
       axes = FALSE, xlab = xlab, ylab = ylab)
plot.new()
make.plot(x[, c(3, 1)]); axis(3); axis(4)
make.plot(x[, c(1, 2)], "Time", "X1"); axis(1); axis(2)
make.plot(x[, c(3, 2)], "X2"); axis(1); axis(4)
title(main = "Van der Pol Oscillator", outer = TRUE)

# Use a dynamic parameter
VdP.sys = '
dxdt[0] = x[1];
dxdt[1] = mu * (1 - x[0] * x[0]) * x[1] - x[0];
' # Vdp.sys
compile_sys("vpol2", VdP.sys, "mu", method = "bsd")
par(mfrow = c(2, 2), mar = rep(1, 4), oma = rep(3, 4), xpd = NA)
for (mu in seq(0.5, 2, len = 4))
{
  vpol2_set_params(mu = mu)
  x = vpol2(rep(1e-4, 2), 100, 0.01)
  make.plot(x[, 2:3]); box()
  title(paste("mu =", round(mu, 2)))
}
title("Van der Pol Oscillator Parameter Sweep", outer = TRUE)
title(xlab = "X1", ylab = "X2", line = 0, outer = TRUE)

# Stiff example from odeint docs
Stiff.sys = '
  dxdt[0] = -101.0 * x[0] - 100.0 * x[1];
  dxdt[1] = x[0];
' # Stiff.sys
cat(JacobianCpp(Stiff.sys))
J(0, 0) = -101;
J(0, 1) = -100;
J(1, 0) = 1;
J(1, 1) = 0;
dfdt[0] = 0.0;
dfdt[1] = 0.0;
compile_implicit("stiff", Stiff.sys)
x = stiff(c(2, 1), 5, 0.001)
plot(x[, 1:2], type = "l", lwd = 2, col = "steelblue")
lines(x[, c(1, 3)], lwd = 2, col = "darkred")

Numerical Integration

pracma() Package

Description: The pracma package offers a wide range of numerical methods, including tools for numerical integration. It is useful for performing one-dimensional integration as well as handling more complex mathematical tasks.

library(pracma)

Attaching package: 'pracma'
The following object is masked from 'package:deSolve':

    rk4
# Define the function to integrate
f <- function(x) sin(x)

# Perform numerical integration
result <- integral(f, 0, pi)
print(result)
[1] 2

cubature() Package

Description: The cubature package is designed for adaptive multivariate integration over hypercubes. It is particularly useful for higher-dimensional integration problems and can handle a variety of complex integrals efficiently.

library(cubature)

# Define the multivariate function to integrate
f <- function(x) x[1]^2 + x[2]^2

# Perform numerical integration over [0, 1] x [0, 1]
result <- adaptIntegrate(f, lowerLimit = c(0, 0), upperLimit = c(1, 1))
print(result$integral)
[1] 0.6666667

Base R integrate Function

Description: The integrate function in base R is a straightforward tool for performing numerical integration of one-dimensional functions. It is easy to use and provides reliable results for a wide range of integrals.

# Define the function to integrate
f <- function(x) x^2

# Perform numerical integration
result <- integrate(f, 0, 1)
print(result$value)
[1] 0.3333333

These notes should provide your readers with a clear understanding of the recommended packages for solving stiff differential equations and performing numerical integration in R, along with example codes and references for further reading.

R has several packages that are well-regarded for integrating and solving differential equations. Here are some of the most popular ones:

  1. deSolve():

    • This package provides functions to solve initial value problems of differential equations, including ordinary differential equations (ODEs), partial differential equations (PDEs), and delay differential equations (DDEs).

    • It supports both stiff and non-stiff systems and offers a variety of solvers.

    • Example usage

      library(deSolve)

# Define the ODE function
model <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    dN <- r * N
    return(list(c(dN)))
  })
}

# Initial state and parameters
state <- c(N = 1)
parameters <- c(r = 0.1)
times <- seq(0, 100, by = 1)

# Solve the ODE
out <- ode(y = state, times = times, func = model, parms = parameters)

  2. PBSddesolve():

  • This package is designed for solving delay differential equations (DDEs) and ordinary differential equations (ODEs). It is useful if your system involves delays.

  • Example usage

library(PBSddesolve)

-----------------------------------------------------------
PBSddesolve 1.13.4 -- Copyright (C) 2007-2024 Fisheries and Oceans Canada
(based on solv95 by Simon Wood)

A complete user guide 'PBSddesolve-UG.pdf' is located at 
C:/Users/maths/AppData/Local/R/win-library/4.3/PBSddesolve/doc/PBSddesolve-UG.pdf

Demos include 'blowflies', 'cooling', 'icecream', and 'lorenz'
They can be run two ways:

1. Using 'utils' package 'demo' function, run command
> demo(icecream, package='PBSddesolve', ask=FALSE)

2. Using package 'PBSmodelling', run commands
> require(PBSmodelling); runDemos()
and choose PBSddesolve and then one of the four demos.

Packaged on 2024-01-04
Pacific Biological Station, Nanaimo

All available PBS packages can be found at
https://github.com/pbs-software
-----------------------------------------------------------
# Define the DDE function
model <- function(t, y, parms) {
  dy <- parms["r"] * y[1]
  return(list(dy))
}

# Initial state and parameters
state <- c(N = 1)
parameters <- c(r = 0.1)
times <- seq(0, 100, by = 1)

# Solve the DDE
out <- dde(y = state, times = times, func = model, parms = parameters)
  1. ReacTran():

  • Useful for solving reaction-transport equations, especially in one, two, or three dimensions. It builds on the deSolve package.

  • Example usage

    library(ReacTran)
    Loading required package: rootSolve
    
Attaching package: 'rootSolve'
    The following objects are masked from 'package:pracma':

    gradient, hessian
    Loading required package: shape
    library(deSolve)

# Example setup for ReacTran

Solving integral equations in R can be accomplished using several packages, each with its specific strengths for different types of integral equations. Here are a few notable ones:

  1. pracma():

    • This package provides numerical tools for a variety of mathematical operations, including solving integral equations.

    • Example usage

      library(pracma)

# Define the integral equation
f <- function(x) x^2

# Solve the integral from 0 to 1
result <- integral(f, 0, 1)
print(result)
      [1] 0.3333333

  2. DEoptim:

  • This package is primarily used for optimization, but it can also handle integral equations through its optimization routines by minimizing the difference between the left and right sides of the integral equation.

  • Example usage:

    library(DEoptim)
    Loading required package: parallel
    
DEoptim package
Differential Evolution algorithm in R
Authors: D. Ardia, K. Mullen, B. Peterson and J. Ulrich
    # Define the integral equation
f <- function(x) x^2
integral_eq <- function(a) abs(integral(f, 0, a) - target_value)

# Use DEoptim to solve for a
target_value <- 1
result <- DEoptim(integral_eq, lower = 0, upper = 2)
    Iteration: 1 bestvalit: 0.040710 bestmemit:    1.461561
Iteration: 2 bestvalit: 0.008380 bestmemit:    1.438209
Iteration: 3 bestvalit: 0.008380 bestmemit:    1.438209
Iteration: 4 bestvalit: 0.008380 bestmemit:    1.438209
Iteration: 5 bestvalit: 0.008380 bestmemit:    1.438209
Iteration: 6 bestvalit: 0.008380 bestmemit:    1.438209
Iteration: 7 bestvalit: 0.008380 bestmemit:    1.438209
Iteration: 8 bestvalit: 0.008380 bestmemit:    1.438209
Iteration: 9 bestvalit: 0.008380 bestmemit:    1.438209
Iteration: 10 bestvalit: 0.001302 bestmemit:    1.442875
Iteration: 11 bestvalit: 0.000350 bestmemit:    1.442081
Iteration: 12 bestvalit: 0.000350 bestmemit:    1.442081
Iteration: 13 bestvalit: 0.000004 bestmemit:    1.442251
Iteration: 14 bestvalit: 0.000004 bestmemit:    1.442251
Iteration: 15 bestvalit: 0.000004 bestmemit:    1.442251
Iteration: 16 bestvalit: 0.000004 bestmemit:    1.442251
Iteration: 17 bestvalit: 0.000004 bestmemit:    1.442251
Iteration: 18 bestvalit: 0.000004 bestmemit:    1.442251
Iteration: 19 bestvalit: 0.000004 bestmemit:    1.442251
Iteration: 20 bestvalit: 0.000004 bestmemit:    1.442251
Iteration: 21 bestvalit: 0.000004 bestmemit:    1.442251
Iteration: 22 bestvalit: 0.000004 bestmemit:    1.442251
Iteration: 23 bestvalit: 0.000002 bestmemit:    1.442251
Iteration: 24 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 25 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 26 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 27 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 28 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 29 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 30 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 31 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 32 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 33 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 34 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 35 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 36 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 37 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 38 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 39 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 40 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 41 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 42 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 43 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 44 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 45 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 46 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 47 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 48 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 49 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 50 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 51 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 52 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 53 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 54 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 55 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 56 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 57 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 58 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 59 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 60 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 61 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 62 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 63 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 64 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 65 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 66 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 67 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 68 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 69 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 70 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 71 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 72 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 73 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 74 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 75 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 76 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 77 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 78 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 79 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 80 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 81 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 82 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 83 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 84 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 85 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 86 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 87 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 88 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 89 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 90 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 91 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 92 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 93 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 94 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 95 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 96 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 97 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 98 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 99 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 100 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 101 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 102 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 103 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 104 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 105 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 106 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 107 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 108 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 109 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 110 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 111 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 112 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 113 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 114 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 115 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 116 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 117 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 118 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 119 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 120 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 121 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 122 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 123 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 124 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 125 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 126 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 127 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 128 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 129 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 130 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 131 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 132 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 133 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 134 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 135 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 136 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 137 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 138 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 139 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 140 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 141 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 142 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 143 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 144 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 145 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 146 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 147 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 148 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 149 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 150 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 151 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 152 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 153 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 154 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 155 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 156 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 157 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 158 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 159 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 160 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 161 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 162 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 163 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 164 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 165 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 166 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 167 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 168 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 169 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 170 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 171 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 172 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 173 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 174 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 175 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 176 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 177 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 178 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 179 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 180 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 181 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 182 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 183 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 184 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 185 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 186 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 187 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 188 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 189 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 190 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 191 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 192 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 193 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 194 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 195 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 196 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 197 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 198 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 199 bestvalit: 0.000000 bestmemit:    1.442250
Iteration: 200 bestvalit: 0.000000 bestmemit:    1.442250
    print(result$optim$bestmem)
       par1 
1.44225
  1. integrate (base R):

  • The integrate function in base R can be used to perform numerical integration, which is essential for solving certain integral equations.

  • Example usage

    f <- function(x) x^2

# Solve the integral from 0 to 1
result <- integrate(f, 0, 1)
print(result$value)
    [1] 0.3333333
  1. Cubature:

  • This package is used for adaptive multivariate integration over hypercubes, which can be useful for higher-dimensional integral equations.

  • Example usage

    library(cubature)

# Define the multivariate integral function
f <- function(x) x[1]^2 + x[2]^2

# Solve the integral over [0, 1] x [0, 1]
result <- adaptIntegrate(f, lowerLimit = c(0, 0), upperLimit = c(1, 1))
print(result$integral)
    [1] 0.6666667
  1. integral (pracma):

  • In addition to the general pracma package, the specific integral function can be quite handy.

  • Example usage

    library(pracma)

# Define the integral equation function
f <- function(x) sin(x)

# Perform the integration over the range 0 to pi
result <- integral(f, 0, pi)
print(result)
    [1] 2

Each of these packages can handle different types of integral equations, and their selection depends on the specific problem and complexity of the equations you are dealing with. For many general purposes, the integrate function in base R and the pracma package’s integral function are very effective and straightforward to use.

Cran Documentation References:

  1. deSolve Package Documentation

  2. odeintr Package Documentation

  3. pracma Package Documentation

  4. cubature Package Documentation

  5. cubature: Adaptive Multivariate Integration over Hypercubes