Merge pull request #779 from hurak/master

SDE Examples: minor rewording of explanation of "diagonal vs non-diagonal noise"

Author: ChrisRackauckas

docs/make.jl

@@ -39,7 +39,9 @@ makedocs(
         "https://www.wolframalpha.com/input/?i=u%27%3D-sqrt%28u%29",
         "https://www.mathworks.com/help/simulink/gui/absolutetolerance.html",
         "https://www.mathworks.com/help/matlab/math/choose-an-ode-solver.html",
-        "https://journals.ametsoc.org/view/journals/atsc/20/2/1520-0469_1963_020_0130_dnf_2_0_co_2.xml"
+        "https://journals.ametsoc.org/view/journals/atsc/20/2/1520-0469_1963_020_0130_dnf_2_0_co_2.xml",
+        "https://mathematica.stackexchange.com/questions/40122/help-to-plot-poincar%C3%A9-section-for-double-pendulum",
+        "https://github.com/SciML/DiffEqProblemLibrary.jl/blob/master/lib/SDEProblemLibrary/src/SDEProblemLibrary.jl",
     ],
     doctest = false, clean = true,
     warnonly = [:missing_docs],

docs/src/tutorials/sde_example.md

@@ -142,49 +142,45 @@ In general, a system of SDEs
 du = f(u,p,t)dt + g(u,p,t)dW,
 ```
 
-where `g` is now a matrix of values, is numerically integrated in the
-same way as ODEs. A common scenario is when we have diagonal noise, which
-is the default for DifferentialEquations.jl. Physically this means that
-every variable in the system gets a different random kick. Consequently, `g` is a
-diagonal matrix and we can handle this in a simple manner by defining
-the deterministic part `f(du,u,p,t)` and the stochastic part
-`g(du2,u,p,t)` as in-place functions.
-
-Consider for example a stochastic variant of the Lorenz equations, where we introduce a
-simple additive noise to each of `x,y,z`, which is simply `3*N(0,dt)`. Here, `N` is the normal
-distribution and `dt` is the time step. This is done as follows:
+where `u` is now a vector of variables, `f` is a vector, and `g` is a matrix, is numerically integrated in the same way as ODEs. A common scenario, which is the default for DifferentialEquations.jl, is that every variable in the system gets a different random kick. This is the case when `g` is a diagonal matrix. Correspondingly, we say that we have a diagonal noise. 
+
+We handle this in a simple manner by defining the deterministic part `f!(du,u,p,t)` and the stochastic part
+`g!(du2,u,p,t)` as in-place functions, but note that our convention is that the function `g!` only defines and modifies the diagonal entries of `g` matrix.
+
+As an example, we consider a stochastic variant of the Lorenz equations, where we add to each of `u₁, u₂, u₃` their own simple noise `3*N(0,dt)`. Here, `N` is the normal distribution and `dt` is the time step. This is done as follows:
 
 ```@example sde2
 using DifferentialEquations
 using Plots
-function lorenz!(du, u, p, t)
+
+function f!(du, u, p, t)
     du[1] = 10.0 * (u[2] - u[1])
     du[2] = u[1] * (28.0 - u[3]) - u[2]
     du[3] = u[1] * u[2] - (8 / 3) * u[3]
 end
 
-function σ_lorenz!(du, u, p, t)
+function g!(du, u, p, t)  # It actually represents a diagonal matrix [3.0 0 0; 0 3.0 0; 0 0 3.0]
     du[1] = 3.0
     du[2] = 3.0
     du[3] = 3.0
 end
 
-prob_sde_lorenz = SDEProblem(lorenz!, σ_lorenz!, [1.0, 0.0, 0.0], (0.0, 10.0))
+prob_sde_lorenz = SDEProblem(f!, g!, [1.0, 0.0, 0.0], (0.0, 10.0))
 sol = solve(prob_sde_lorenz)
 plot(sol, idxs = (1, 2, 3))
 ```
 
 Note that it's okay for the noise function to mix terms. For example
 
 ```@example sde2
-function σ_lorenz!(du, u, p, t)
+function g!(du, u, p, t)
     du[1] = sin(u[3]) * 3.0
     du[2] = u[2] * u[1] * 3.0
     du[3] = 3.0
 end
 ```
 
-is a valid noise function, which will once again give diagonal noise by `du2.*W`.
+is a valid noise function.
 
 ## Example 3: Systems of SDEs with Scalar Noise
 
@@ -200,12 +196,12 @@ let's solve a linear SDE with scalar noise using a high order algorithm:
 ```@example sde3
 using DifferentialEquations
 using Plots
-f(du, u, p, t) = (du .= u)
-g(du, u, p, t) = (du .= u)
+f!(du, u, p, t) = (du .= u)
+g!(du, u, p, t) = (du .= u)
 u0 = rand(4, 2)
 
 W = WienerProcess(0.0, 0.0, 0.0)
-prob = SDEProblem(f, g, u0, (0.0, 1.0), noise = W)
+prob = SDEProblem(f!, g!, u0, (0.0, 1.0), noise = W)
 sol = solve(prob, SRIW1())
 plot(sol)
 ```
@@ -216,7 +212,7 @@ In the previous examples we had diagonal noise, that is a vector of random numbe
 `dW` whose size matches the output of `g` where the noise is applied element-wise,
 and scalar noise where a single random variable is applied to all dependent variables.
 However, a more general type of noise allows for the terms to linearly mixed via `g`
-being a matrix.
+being a general nondiagonal matrix.
 
 Note that nonlinear mixings are not SDEs but fall under the more general class of
 random ordinary differential equations (RODEs) which have a
@@ -228,8 +224,8 @@ is `g(u,p,t)*dW`, the matrix multiplication. For example, we can do the followin
 
 ```@example sde4
 using DifferentialEquations
-f(du, u, p, t) = du .= 1.01u
-function g(du, u, p, t)
+f!(du, u, p, t) = du .= 1.01u
+function g!(du, u, p, t)
     du[1, 1] = 0.3u[1]
     du[1, 2] = 0.6u[1]
     du[1, 3] = 0.9u[1]
@@ -239,10 +235,10 @@ function g(du, u, p, t)
     du[2, 3] = 0.3u[2]
     du[2, 4] = 1.8u[2]
 end
-prob = SDEProblem(f, g, ones(2), (0.0, 1.0), noise_rate_prototype = zeros(2, 4))
+prob = SDEProblem(f!, g!, ones(2), (0.0, 1.0), noise_rate_prototype = zeros(2, 4))
 ```
 
-In our `g` we define the functions for computing the values of the matrix.
+In our `g!` we define the functions for computing the values of the matrix.
 We can now think of the SDE that this solves as the system of equations
 
 ```math
@@ -259,7 +255,7 @@ the same random number in the first and second SDEs.
     noise. This is discussed [in the SDE solvers page](@ref sde_solve).
 
 The matrix itself is determined by the keyword argument `noise_rate_prototype` in the `SDEProblem`
-constructor. This is a prototype for the type that `du` will be in `g`. This can
+constructor. This is a prototype for the type that `du` will be in `g!`. This can
 be any `AbstractMatrix` type. Thus, we can define the problem as
 
 ```@example sde4
@@ -271,18 +267,18 @@ A[1, 4] = 1
 A[2, 4] = 1
 A = sparse(A)
 
-# Make `g` write the sparse matrix values
-function g(du, u, p, t)
+# Make `g!` write the sparse matrix values
+function g!(du, u, p, t)
     du[1, 1] = 0.3u[1]
     du[1, 4] = 0.12u[2]
     du[2, 4] = 1.8u[2]
 end
 
-# Make `g` use the sparse matrix
-prob = SDEProblem(f, g, ones(2), (0.0, 1.0), noise_rate_prototype = A)
+# Make `g!` use the sparse matrix
+prob = SDEProblem(f!, g!, ones(2), (0.0, 1.0), noise_rate_prototype = A)
 ```
 
-and now `g(u,p,t)` writes into a sparse matrix, and `g(u,p,t)*dW` is sparse matrix
+and now `g!(u,p,t)` writes into a sparse matrix, and `g!(u,p,t)*dW` is sparse matrix
 multiplication.
 
 ## Example 4: Colored Noise
@@ -292,7 +288,7 @@ In that portion of the docs, it is shown how to define your own noise process
 `my_noise`, which can be passed to the SDEProblem
 
 ```julia
-SDEProblem(f, g, u0, tspan, noise = my_noise)
+SDEProblem(f!, g!, u0, tspan, noise = my_noise)
 ```
 
 Note that general colored noise problems are only compatible with the `EM` and `EulerHeun` methods.
@@ -311,11 +307,11 @@ dW_1 dW_2 = ρ dt
 In this problem, we have a diagonal noise problem given by:
 
 ```@example sde4
-function f(du, u, p, t)
+function f!(du, u, p, t)
     du[1] = μ * u[1]
     du[2] = κ * (Θ - u[2])
 end
-function g(du, u, p, t)
+function g!(du, u, p, t)
     du[1] = √u[2] * u[1]
     du[2] = σ * √u[2]
 end
@@ -339,7 +335,7 @@ heston_noise = CorrelatedWienerProcess!(Γ, tspan[1], zeros(2), zeros(2))
 This is then used to build the SDE:
 
 ```@example sde4
-SDEProblem(f, g, ones(2), tspan, noise = heston_noise)
+SDEProblem(f!, g!, ones(2), tspan, noise = heston_noise)
 ```
 
 Of course, to fully define this problem, we need to define our constants. Constructors