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

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.

• 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. Commented Jul 18, 2019 at 15:29
• 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. Commented Jul 18, 2019 at 15:42
• Any of the books in the answers to this question will have the information you need: math.stackexchange.com/questions/512432/… Commented Jul 18, 2019 at 17:44
• 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 Commented Jul 19, 2019 at 6:17
• 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. Commented Jul 20, 2019 at 18:27

## 1 Answer

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.

• 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? Commented Jul 20, 2019 at 18:30
• 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. Commented Jul 21, 2019 at 0:41
• 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. Commented Jul 21, 2019 at 17:57
• 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? Commented Jul 21, 2019 at 18:49
• 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. Commented Jul 21, 2019 at 22:17