How to use Newton-Raphson method to handle nonlinear terms in coupled system of PDEs?

Question

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.

Answer 1

The Nonlinear Schroedinger Equation is special. Instead of applying a Newton-Raphson method, it is easier to use an operator splitting scheme that uses the particular form of the nonlinearity because the nonlinear part of the equation allows for an analytical solution.

I started writing a tutorial program for the deal.II library a while back that illustrates this. I never managed to finish it, but most of the documentation is actually ready: https://github.com/bangerth/dealii/blob/step-58/examples/step-58/doc/intro.dox I suspect it might be difficult to read this document at the location of the link. Contact me offline if you want a better form.