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?


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.

  • $\begingroup$ The topic you're asking about is numerical solution of nonlinear equations, and is a well-established field of research. Newton-Rhapson is the most basic method, for scalar equations. If you're looking for software that will solve your equations, let us know what language you are using. If you're looking for understanding, any introductory text on numerical analysis will have a chapter on this topic. $\endgroup$ – David Ketcheson Jul 18 '19 at 15:29
  • $\begingroup$ I'm looking for understanding. If you could suggest me some literature I'd ve very much appreciated. I've checked some, but my problem is to handle the system ( the two vectors of variables /rho and /psi). This is not very clear... What i mean, specifically, is that i don't know how to apply newton's method to a system. For example, should i start by saying /rho=/rho_0+/delta_1 and /psi=/psi_0+/delta_2 and obtain a linear system for the deltas? How would i compute the jacobian in this case? I'd be very much appreciated if someone could go through the method in my specific case. $\endgroup$ – Rui Martins Jul 18 '19 at 15:42
  • $\begingroup$ Any of the books in the answers to this question will have the information you need: math.stackexchange.com/questions/512432/… $\endgroup$ – David Ketcheson Jul 18 '19 at 17:44
  • $\begingroup$ How much of this do you have to code yourself? There are very robust methods for solving large, nonlinear, and stiff systems, e.g. DifferentialEquations.jl. These incorporate Newton's Method in very efficient ways and handle step sizes and mass matrices like in your problem. Even if you can't use them, it may be worth checking out their algorithms and seeing what you can learn $\endgroup$ – whpowell96 Jul 19 '19 at 6:17
  • $\begingroup$ Thanks for all the help! In the last couple of days I've been trying to read the literature suggested and trying to implement the algorithm in Matlab by myself . The idea was to code as much as possible myself, so I could learn about FEM and Newton-Rapshon and at the same time solve this equation, that is the subject of my thesis. I'll make another post with what I've been doing soon, so I'll get some help from you guys. $\endgroup$ – Rui Martins Jul 20 '19 at 18:27

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.

| cite | improve this answer | |
  • $\begingroup$ Hello! Yeah, I knew about this scheme and I've been solving the NLSE with that up until now. However, I'd like to not only learn about FEM, but later on to implement mesh adaptivity. For that operator splitting scheme you'd need a regular mesh if I'm not mistaken. Would you say that you can't handle the non-linear term in the NLSE with Newton-Raphson? $\endgroup$ – Rui Martins Jul 20 '19 at 18:30
  • $\begingroup$ Of course you can handle the nonlinearity with a Newton-Raphson scheme. It's just unnecessary work. As for mesh adaptivity: I see no reason why it can't be used together with the operator splitting scheme. The operator splitting is simply a way to solve a single time step. What you do between time steps (e.g., adapting the mesh) remains unaffected. $\endgroup$ – Wolfgang Bangerth Jul 21 '19 at 0:41
  • $\begingroup$ Yeah, but the operator splitting scheme involves DFT, and I'm not quite certain, but you do need a regular mesh to perform them right? The idea was to have an irregular spatial mesh focusing on the difficult parts. Nevertheless I'd like to be able to solve the equation with FE, without worrying about efficiency. $\endgroup$ – Rui Martins Jul 21 '19 at 17:57
  • $\begingroup$ I'm new into numerical methods. Maybe what i just said, thinking a little bit about it, does not make too much sense. So you say that mesh adaptivity, even by making it very irregular ( coarser at some regions, well defined at others), is applicable with the operator splitting scheme? $\endgroup$ – Rui Martins Jul 21 '19 at 18:49
  • $\begingroup$ The operator splitting just requires you to solve two sub-steps. One of those (the nonlinear one) is a local operator on each node. The other one requires you to solve a PDE. It's true that you can do that using the DFT if you want -- and indeed, in that case you probably want a uniform mesh. But you can also solve that PDE like we do for all other equations, by building a linear system and then solving it. $\endgroup$ – Wolfgang Bangerth Jul 21 '19 at 22:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.