| Title: | Monte Carlo Solution of First Order Differential Equations |
|---|---|
| Description: | Two functions for simulating the solution of initial value problems of the form g'(x) = G(x, g) with g(x0) = g0. One is an acceptance-rejection method. The other is a method based on the Mean Value Theorem. |
| Authors: | W.J. Braun [aut, cre] |
| Maintainer: | W.J. Braun <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 1.1 |
| Built: | 2024-09-20 10:07:32 UTC |
| Source: | CRAN |
Monte Carlo solution of first order differential equations and confidence intervals for the solutions.
ODE.MVT can solve any first order differential equation of the form g' = F(x, g) with initial condition g(x0) = g0, provided conditions are satisfied for existence and uniqueness of the solution. Confidence intervals for the solution can be obtained through use of ODE.MVT.agg. An alternate method due to Akhtar et al is implemented in ODE.ADA.
Author: W.J. Braun
Maintainer: W.J. Braun
Akhtar, M. N., Durad, M. H., and Ahmed, A. (2015). Solving initial value ordinary differential equations by Monte Carlo method. Proc. IAM, 4:149-174.
Braun, W. J. (2024) Monte Carlo integration of first order ordinary differential equations. Preprint.
deSolve
# Solve g' = F(x, g) on (0, 1] with g(0) = -1/1000001 G <- function(x, g) { -1000*g + 3000 - 2000*exp(-x) } T <- 1 x0 <- 0 g0 <- -1/1000001 nMVT <- 5000 ghat <- ODE.MVT(G, initvalue = g0, endpoint = T, initpoint = x0, Niter = 2, npoints = nMVT) plot(ghat, type = "l")# Solve g' = F(x, g) on (0, 1] with g(0) = -1/1000001 G <- function(x, g) { -1000*g + 3000 - 2000*exp(-x) } T <- 1 x0 <- 0 g0 <- -1/1000001 nMVT <- 5000 ghat <- ODE.MVT(G, initvalue = g0, endpoint = T, initpoint = x0, Niter = 2, npoints = nMVT) plot(ghat, type = "l")
Given g' = G(x, g) and g(x0) = g0 satisfying conditions for existence and uniqueness of solution, a Monte Carlo approximation to the solution is found. The method is essentially a rejection method.
ODE.ADA(G, initvalue, endpoint, X0 = 0, npoints = 1000, M, R, N = 1e+05)ODE.ADA(G, initvalue, endpoint, X0 = 0, npoints = 1000, M, R, N = 1e+05)
G |
a function having two numeric vector arguments. |
initvalue |
a numeric initial value |
endpoint |
a numeric interval endpoint value. |
X0 |
a numeric interval starting point value. |
npoints |
an integer value specifying the number of subintervals to build the approximation on. |
M |
a numeric value specifying an upper bound for F(x, g). |
R |
a numeric value specifying a lower bound for F(x, g). |
N |
an integer value specifying the number of Monte Carlo simulations used for each subinterval. |
A list consisting of
x |
a vector of length npoints, consisting of the x-coordinate of the solution. |
y |
a vector of length npoints, consisting of the y-coordinate of the solution, i.e. g(x). |
Akhtar, M. N., Durad, M. H., and Ahmed, A. (2015). Solving initial value ordinary differential equations by Monte Carlo method. Proc. IAM, 4:149-174.
G <- function(x, g) { -1000*g + sin(x) } gTrue <- function(x) (1000*sin(x) - cos(x))/1000001 xF <- 1; x0 <- 0 g0 <- -1/1000001 # initial value NAkhtar <- 800 nAkhtar <- 2500 MAkhtar <- 1e-3 RAkhtar <- 1e-6 gAkhtar <- ODE.ADA(G, initvalue = g0, endpoint = xF, X0 = x0, npoints = nAkhtar, M=MAkhtar, R = RAkhtar, N = NAkhtar) plot(gAkhtar, type = "l")G <- function(x, g) { -1000*g + sin(x) } gTrue <- function(x) (1000*sin(x) - cos(x))/1000001 xF <- 1; x0 <- 0 g0 <- -1/1000001 # initial value NAkhtar <- 800 nAkhtar <- 2500 MAkhtar <- 1e-3 RAkhtar <- 1e-6 gAkhtar <- ODE.ADA(G, initvalue = g0, endpoint = xF, X0 = x0, npoints = nAkhtar, M=MAkhtar, R = RAkhtar, N = NAkhtar) plot(gAkhtar, type = "l")
Given g' = G(x, g) and g(x0) = g0 satisfying conditions for existence and uniqueness of solution, a Monte Carlo approximation to the solution is found. The method is essentially a Monte Carlo version of Euler and Runge-Kutta type methods.
ODE.MVT(G, initvalue, endpoint, initpoint = 0, Niter = 2, npoints = 1000)ODE.MVT(G, initvalue, endpoint, initpoint = 0, Niter = 2, npoints = 1000)
G |
a function having two arguments. |
initvalue |
a numeric initial value |
endpoint |
a numeric interval endpoint value. |
initpoint |
a numeric interval starting point value. |
Niter |
an integer value specifying the number of iterations at each step. |
npoints |
an integer value specifying the number of subintervals to build the approximation on. |
A list consisting of
x |
a vector of length npoints, consisting of the x-coordinate of the solution. |
y |
a vector of length npoints, consisting of the y-coordinate of the solution, i.e. g(x). |
Braun, W.J.
Braun, W.J. (2024) Preprint.
# Initial Value Problem: G <- function(x, g) { -1000*g + sin(x) } g0 <- -1/1000001; x0 <- 0 # initial condition xF <- 1 # interval endpoint nMVT <- 1000 # number of subintervals used # Monte Carlo solution based on one iteration: ghat1 <- ODE.MVT(G, initvalue = g0, endpoint = xF, initpoint = x0, Niter = 1, npoints = nMVT) # Monte Carlo solution based on five iterations: ghat5 <- ODE.MVT(G, initvalue = g0, endpoint = xF, initpoint = x0, Niter = 5, npoints = nMVT) gTrue <- function(x) (1000*sin(x) - cos(x))/1000001 # true solution oldpar <- par(mfrow=c(1,3)) curve(gTrue(x)) # lines(ghat1, col = 2, lty = 2, ylab="g(x)") plot(abs(gTrue(ghat1$x) - ghat1$y) ~ ghat1$x, xlab = "x", ylab = "Absolute Error", type = "l", ylim = c(0, 2e-6), main="1 Iteration") plot(abs(gTrue(ghat5$x) - ghat5$y) ~ ghat5$x, xlab = "x", ylab = "Absolute Error", type = "l", ylim = c(0, 2e-6), main="5 Iterations") par(oldpar)# Initial Value Problem: G <- function(x, g) { -1000*g + sin(x) } g0 <- -1/1000001; x0 <- 0 # initial condition xF <- 1 # interval endpoint nMVT <- 1000 # number of subintervals used # Monte Carlo solution based on one iteration: ghat1 <- ODE.MVT(G, initvalue = g0, endpoint = xF, initpoint = x0, Niter = 1, npoints = nMVT) # Monte Carlo solution based on five iterations: ghat5 <- ODE.MVT(G, initvalue = g0, endpoint = xF, initpoint = x0, Niter = 5, npoints = nMVT) gTrue <- function(x) (1000*sin(x) - cos(x))/1000001 # true solution oldpar <- par(mfrow=c(1,3)) curve(gTrue(x)) # lines(ghat1, col = 2, lty = 2, ylab="g(x)") plot(abs(gTrue(ghat1$x) - ghat1$y) ~ ghat1$x, xlab = "x", ylab = "Absolute Error", type = "l", ylim = c(0, 2e-6), main="1 Iteration") plot(abs(gTrue(ghat5$x) - ghat5$y) ~ ghat5$x, xlab = "x", ylab = "Absolute Error", type = "l", ylim = c(0, 2e-6), main="5 Iterations") par(oldpar)
Given g' = G(x, g) and g(x0) = g0 satisfying conditions for existence and uniqueness
of solution, a Monte Carlo approximation to the solution is found. The method is
essentially a Monte Carlo version of Euler and Runge-Kutta type methods.
Repeated calls are made to the ODE.MVT function to obtain replicated estimates
of the solution.
ODE.MVT.agg(G, initvalue, endpoint, initpoint = 0, Niter = 2, npoints = 1000, nsims = 10)ODE.MVT.agg(G, initvalue, endpoint, initpoint = 0, Niter = 2, npoints = 1000, nsims = 10)
G |
a function having two numeric vector arguments. |
initvalue |
a numeric initial value |
endpoint |
a numeric interval endpoint value. |
initpoint |
a numeric interval starting point value. |
Niter |
an integer value specifying the number of iterations at each step. |
npoints |
an integer value specifying the number of subintervals to build the approximation on. |
nsims |
an integer value specifying the number of replicated solutions required. |
A list consisting of
x |
a vector of length npoints, consisting of the x-coordinate of the solution. |
y |
a vector of length npoints, consisting of the average of replicated y-coordinates of the solution, i.e. g(x). |
stdev |
a vector of length npoints, consisting of the standard deviation estimates of the solution estimate at each point. |
W.J. Braun
G <- function(x, g) { -1000*g + sin(x) } g0 <- -1/1000001; x0 <- 0 # initial condition xF <- 1 # interval endpoint nMVT <- 1000 # number of subintervals used # 10 reps of Monte Carlo solution based on two iterations: ghat2 <- ODE.MVT.agg(G, initvalue = g0, endpoint = xF, initpoint = x0, Niter = 2, npoints = nMVT, nsims = 10) # 10 reps of Monte Carlo solution based on three iterations: ghat3 <- ODE.MVT.agg(G, initvalue = g0, endpoint = xF, initpoint = x0, Niter = 3, npoints = nMVT, nsims = 10) gTrue <- function(x) (1000*sin(x) - cos(x))/1000001 # true solution oldpar <- par(mfrow=c(1,2)) curve(gTrue(x), from = 0.3, to = 0.3015) # graph of solution on small interval lines(ghat2, col = 4, lty = 4, ylab="g(x)", lwd = 2) # estimated MC solution - 2 iterations lines(ghat2$y - 1.96*ghat2$stdev ~ ghat2$x, lty = 3, col = 4) # lower conf. limits lines(ghat2$y + 1.96*ghat2$stdev ~ ghat2$x, lty = 3, col = 4) # upper conf. limits curve(gTrue(x), from = 0.3, to = 0.3015) # graph of solution on small interval lines(ghat3, col = 4, lty = 4, ylab="g(x)", lwd = 2) # estimated MC solution - 3 iterations lines(ghat3$y - 1.96*ghat3$stdev ~ ghat3$x, lty = 3, col = 4) # lower conf. limits lines(ghat3$y + 1.96*ghat3$stdev ~ ghat3$x, lty = 3, col = 4) # upper conf. limits par(oldpar)G <- function(x, g) { -1000*g + sin(x) } g0 <- -1/1000001; x0 <- 0 # initial condition xF <- 1 # interval endpoint nMVT <- 1000 # number of subintervals used # 10 reps of Monte Carlo solution based on two iterations: ghat2 <- ODE.MVT.agg(G, initvalue = g0, endpoint = xF, initpoint = x0, Niter = 2, npoints = nMVT, nsims = 10) # 10 reps of Monte Carlo solution based on three iterations: ghat3 <- ODE.MVT.agg(G, initvalue = g0, endpoint = xF, initpoint = x0, Niter = 3, npoints = nMVT, nsims = 10) gTrue <- function(x) (1000*sin(x) - cos(x))/1000001 # true solution oldpar <- par(mfrow=c(1,2)) curve(gTrue(x), from = 0.3, to = 0.3015) # graph of solution on small interval lines(ghat2, col = 4, lty = 4, ylab="g(x)", lwd = 2) # estimated MC solution - 2 iterations lines(ghat2$y - 1.96*ghat2$stdev ~ ghat2$x, lty = 3, col = 4) # lower conf. limits lines(ghat2$y + 1.96*ghat2$stdev ~ ghat2$x, lty = 3, col = 4) # upper conf. limits curve(gTrue(x), from = 0.3, to = 0.3015) # graph of solution on small interval lines(ghat3, col = 4, lty = 4, ylab="g(x)", lwd = 2) # estimated MC solution - 3 iterations lines(ghat3$y - 1.96*ghat3$stdev ~ ghat3$x, lty = 3, col = 4) # lower conf. limits lines(ghat3$y + 1.96*ghat3$stdev ~ ghat3$x, lty = 3, col = 4) # upper conf. limits par(oldpar)