# 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. – David Ketcheson Jul 18 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. – Rui Martins Jul 18 at 15:42
• Any of the books in the answers to this question will have the information you need: math.stackexchange.com/questions/512432/… – David Ketcheson Jul 18 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 – whpowell96 Jul 19 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. – Rui Martins Jul 20 at 18:27