Hello!
Let’s say I want to solve numerically the following stochastic partial differential equation:
\begin{equation}
du(x,y,t)=\nabla^{2}u(x,y,t) dt+k\cdot dW^{r}(x,y,t)+i k\cdot dW^{c}(x,y,t))
\end{equation},
where u is a complex function in 2d, i is the imaginary unit, and k is a constant.
As one can see, we have a complex noise. In particular, the Winer process has real part W^{r} and imaginary part W^{c}, which are independent random variables.
Now, in order to solve the equation, I discretize in space, such that I obtain the following system of equations:
du[i,j]=F_{i,j}*u[i,j]dt+k*(dW^{r}[i,j]+idW^{c}[i,j])
I didn’t explicitely write the drift part (the laplacian) to make it shorter, and also because I want to concentrate on the stochastic part, as I think it is the cause of my problems.
So let me rewrite just the stochastic part of the system of equations:
du[i,j]=k\cdot dWr[i,j]+i k \cdot dWc[i,j]
(the “i” multiplying k outside the squared brackets is the imaginary unit and not an index!).
Now, reading the page of DifferentialEquations.jl ()DifferentialEquations.jl documentation, I understood that julia handles the stochastic part of an SDE in the following form:
\begin{equation}
du_{stochastic}=G*dW
\end{equation}
where, if u is a vector, then G is a matrix such that the matrix-vector product G*dW gives the stochastic change of the solution du. (So I guess the solver multiplies G by an array of independent noises dW)
Now, in my case u is a matrix, so in order to put myself in the same scenario described in the documentation, I reshaped u from a matrix (N \times M) into a column vector (N M \times 1) ( usually I have N=M=300).
Having done that, the G matrix in my case will be N M \times 2N M, such that, if I multiply it by a vector of independent noises dW of dimension 2N M \times 1, I get a vector of dimension N M \times 1 corresponding to du. In particular, by looking at my system of equations, my G matrix will have the following form:
\begin{equation}
g=\begin{bmatrix}
k & i \cdot k & 0 & 0 & 0 & 0 & \dots \\
0 & 0 & k & i \cdot k& 0 & 0 & \dots \\
0 & 0 & 0 & 0 & k & i \cdot k & \dots \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \\
\end{bmatrix}
\end{equation}.
Therefore, I implemented the noise prototype as follows:
function create_matrix_G_prototype(imax,jmax,p)
g_matrix_prototype=zeros(Complex{Float64},N*M,2*N*M)
for i in 1:imax*jmax
g_matrix_prototype[i,(2*i)-1]=k
g_matrix_prototype[i,(2*i)]=im*k
end
return sparse(g_matrix_prototype)
end
G_matrix_prototype=create_matrix_G_prototype(imax,jmax,p)
In this way, I can implement the G function to pass to the solver as an in-place function:
function G!(du, u, p, t)
du.=G_matrix_prototype
ende
Then I define the problem and solve it with:
u_0=zeros(Complex{Float64},N*M)
tspan=(0.0,1.0)
problem_r=SDEProblem(F!,G!,u_0,tspan,p,noise_rate_prototype=G_matrix_prototype)
sol_r=solve(problem_r)e
The issue here is that, first of all, the code seems to be terribly slow. Moreover, I get a warning telling me that the dt chosen by the solver is smaller than the dt_{min}:
┌ Warning: dt(8.881784197001252e-16) <= dtmin(8.881784197001252e-16) at t=0.37906220534331775, and step error estimate = 1.1208004104491986. Aborting. There is either an error in your model specification or the true solution is unstable.
This is surprising to me since I didn’t expect the equation to be stiff (but at this point I’m starting to guess it is). Has anyone any idea of what is going wrong? what shall I do about the warning? maybe choose a specific solver for stiff equations?
Also, is there a better approach that I could take for this specific case in order to increase the performance and speed of the code? Is it a good idea to separate real and imaginary parts instead of dealing with complex values?
Thanks!