I'm trying to solve the Nonlinear Schrodinger's Equation (NLSE) in 2D using Finite Elements, but I don't know how to handle the nonlinear term. I suppose I have to apply the Newton-Raphson algortihm to my discretized system of PDEs, but i'm not sure how to proceed.
Note: I just started studying Finite Elements two weeks ago, so I'd appreaciate any advice on how to tackle the problem more efficiently!
The NLSE is,
$i\frac{\partial{u}}{dt}=-\frac{1}{2}\nabla^2u-|u|^2u$, where $u$ is complex valued.
Then $u=r+is$, with $r,s$ real valued functions, and plugging this into the NLSE I obtain a system of coupled PDEs.
$ \left\{ \begin{array}{ll} \partial_tr+\frac{1}{2}\nabla^2s=-(r^2+s^2)s\\ \partial_ts-\frac{1}{2}\nabla^2r=(r^2+s^2)r\\ \end{array} \right. $
Using FE to perform the spatial discretization (and assuming null Dirichlet BC for both the functions and their gradients) , the expansion in hat-functions corresponding to each function $r$ and $s$ is $r_h=\sum_j\rho_j(t)\phi_j$ and $s_h=\sum_j\psi_j(t)\phi_j$, and we get
$ \left\{ \begin{array}{ll} M\dot{\rho}-\frac{1}{2}A\psi=-b_1(\rho,\psi)\\ M\dot{\psi}+\frac{1}{2}A\rho=b_2(\rho,\psi)\\ \end{array} \right. $
where $M$ is the mass matrix, $A$ the stifness matrix and $(b_1)_i=\int_\Omega (r_h^2+s_h^2)s_h \phi_i$, $(b_2)_i=\int_\Omega (r_h^2+s_h^2)r_h \phi_i$.
I've read something sugesting considering the term $|u|^2$ known for every time iteration, and given by $|u_{t-1}|^2$, but after implementing the norm $|u|$ kept growing so I suppose this is not the correct way this non-linearity.
So my question is, how would I apply Newton-Raphson's method to the time discretized version (using backwards Euler) of the above equation, namely,
$ \left\{ \begin{array}{ll} M\psi_t+-M\psi_{t-1}+\frac{1}{2}dtA\rho_t=b_2(\rho_l,\psi_l)dt\\ M\rho_t+-M\rho_{t-1}-\frac{1}{2}dtA\psi_t=-b_1(\rho_l,\psi_l)dt\\ \end{array} \right. $
to be able to handle the nonlinearity?
EDIT
I need some help, because I'm not even sure if I'm on the right track.
Last time, I ended up with the following system of equations (discretized spatially and in time),
$ \left\{ \begin{array}{ll} M\psi_t+\frac{1}{2}dtA\rho_t=b_1(\rho_l,\psi_l)dt+M\psi_{t-1}\\ M\rho_t+\frac{1}{2}dtA\psi_t=b_2(\rho_l,\psi_l)dt+M\rho_{t-1}\\ \end{array} \right. $
NOTE: I redefined $b_1$ and $b_2$, just to be in accordance with my notes, they are $b_{1i}=\int(r_h^2+s_h^2)r_h\phi_i$, $b_{2i}=\int-(r_h^2+s_h^2)s_h\phi_i$.
Now I define $\xi_t=\left[\begin{array}{ll}\rho_t\\\psi_t\end{array}\right]$, and the system becomes,
$\begin{bmatrix} M & \frac{1}{2}dtA\\ -\frac{1}{2}dtA & M \end{bmatrix}\xi_t=dt\left[\begin{array}{ll}b_2\\b_1\end{array}\right]+ \begin{bmatrix} M & 0\\ 0 & M \end{bmatrix}\xi_{t-1}\rightarrow\eta(\xi_t)=0 $ where $\eta(\xi_t)$ is the residual. To apply Newton's Method, I compute the Jacobian $J_{mn}=\frac{\partial\eta(\xi_t)_m}{\partial(\xi_t)_n}=\mathcal{M}_{mn}-dt\frac{\partial\bar{b}_m}{\partial(\xi_t)_n}$, where $\mathcal{M}$ is the left side big matrix in the equation above and $\bar{b}$ is the right side big vector.
By expanding $\bar{b}$ in the basis functions (which may be $b_2$ or $b_1$ depending on the index $m$) I was able to compute $J_{mn}$. For example, lets suppose we're looking at $m=1...n_{nodes}$, $\bar{b}_m=-\sum_{ijk}(\rho_i\rho_j+\psi_i\psi_j)\psi_k\int_{domain}\phi_i\phi_j\phi_k\phi_m$, and so,
$ \frac{\partial\bar{b}_m}{\partial(\xi_t)_n} = -\left\{ \begin{array}{ll} 2\sum_{ik}\rho_i\psi_k\mathcal{I}_{inkm}, \text{for } n=1...n_{nodes} \\ 2\sum_{ij}(\rho_i\rho_j+\psi_i\psi_j)\mathcal{I}_{ijnm}, \text{for } n=n_{nodes}+1...2n_{nodes}\\ \end{array} \right. $,
where $\mathcal{I}_{ijnm}=\int_{domain}\phi_i\phi_j\phi_n\phi_m$.
So, is this correct up until now? And if so, how am I supposed to handle the tensor $\mathcal{I}_{ijnm}$? One thing to notice is that I can't really divide the integral over the domain by a sum of integrals over the elements, since $\mathcal{I}_{ijnm}\neq 0$ when not only all the hat functions $\{i,j,n,m\}$ are inside one element , but also when only one is outside the element.
