Stochastic Finite Element Method

# Stochastic Finite Element Method

This tutorial will introduce you to the functionality for solving SPDEs. Other introductions can be found by checking out DiffEqTutorials.jl. This tutorial assumes you have read the Ordinary Differential Equations tutorial.

For most PDE problem types, we can additionally specify them as a stochastic problem by giving the appropriate optional arguments to the constructor. These arguments are a function σ which is the function multiplied to the Brownian increments \$dW\$, and stochastic, a boolean which we put as true for when the equation is stochastic. Another keyword that is optional is noisetype which specifies the type of noise (the "color" of the noise). By default this is Gaussian (Space-time) White Noise.

The following examples show how to change the tutorial problems into stochastic problems.

## Finite Element Stochastic Poisson Equation

We can solve the same PDE as in the Poisson Tutorial except as the stochastic PDE, \$-Δu=f+gdW\$, with additive space-time white noise by specifying the problem as:

``````f(x) = sin(2π.*x[:,1]).*cos(2π.*x[:,2])
dx = 1//2^(5)
mesh = notime_squaremesh([0 1 0 1],dx,:dirichlet)
prob = PoissonProblem(f,mesh,σ=σ)
sol = solve(prob)``````

This gives the following plot (with adding the deterministic solution from the previous example):

## Finite Element Stochastic Heat Equation

This will solve a nonlinear stochastic heat equation \$u_t=Δu+f+gdW\$ with forcing function `f(u)=.5-u`, noise function `g(u)=100u^2` and initial condition `u0=0`. We would expect this system to rise towards the deterministic steady state `u=2` (but stay in mean a bit below it due to 1st order "Milstein" effects), gaining more noise as it increases. This is specified as follows:

``````f(t,x,u)  = ones(size(x,1)) - .5u
u0_func(x) = zeros(size(x,1))
σ(t,x,u) = 1u.^2
tspan = (0.0,5.0)
dx = 1//2^(3)
dt = 1//2^(11)
mesh = parabolic_squaremesh([0 1 0 1],dx,dt,tspan,:neumann)
u0 = u0_func(mesh.node)
prob = HeatProblem(u0,f,mesh,σ=σ)``````

We use the following code to create an animation of the solution:

``````sol = solve(prob,FEMDiffEqHeatEuler(),save_everystep=true,solver=:LU)
using Plots
animate(sol;zlim=(0,3),cbar=false)``````