# Integrators for Nonlinear/Stiff PDE

It was suggested I ask this question in this section. Anyway:

I have a particular nonlinear PDE of the form

$$u_t(x,t)=iu_{xx}(x,t)+f(x,u(x,t)) \tag{1}$$

Where f is some nonlinear function. With boundary/initial conditions

$$u(0,t)=u(L,t)=0, u(x,0)\approx0 \tag{2}$$

Where $$L$$ is the length of the region considered. The equation is homogeneous so the initial conditions cannot be zero everywhere. Thus the interior points need only be close to zero, just not exactly zero.

My first step to use finite difference on the derivative

$$u_{xx}(x_l)=\frac{1}{\Delta x^2}(u_{l+1}-2u_l+u_{l-1}) \tag{3}$$

Giving now a system of stiff nonlinear ODE

$$\textbf u_t=iD\textbf u+\textbf f(x_l,\textbf u) \tag{4}$$

Where $$D$$ is is a matrix with $$-2$$ on the diagonal, $$1$$ on the first off diagonals, zeros everywhere else, and all divided by the spatial step size squared.

My most successful attempt is to use partial backwards Euler where the derivative is taken at the forward time step and the nonlinear part is taken at the "old" time step, that is:

$$\frac{1}{\Delta t} (\textbf u^{n+1}-\textbf u^{n})=iD\textbf u^{n+1}+\textbf f(x_l,\textbf u^n) \tag{5}$$

Which yields the solution

$$\textbf u^{n+1}=(\frac{1}{\Delta t}I-iD)^{-1}(\frac{1}{\Delta t}\textbf u^n+\textbf f(x_l,\textbf u^n)) \tag{6}$$

Where $$I$$ is the identity matrix. The problem with this method is the solutions distort if I change the time step. I have no idea why.

I tried solving the fully implicit problem via Newton-Raphson but the solutions always converge to zero. I do not know why. It's possible I implemented it incorrectly.

I decided to move on to the method of lines where I would use some kind of established time integrator.

Because of my issues with iterative methods I've looked into non-iterative integrators. The first one I attempted was a Rosebrock method. This, however, gives the incorrect solutions (the Rosebrock solutions are symmetric, the real solution should be asymmetric). And like in my semi-implicit method, changing the time step also changes the solutions. Next, I attempted exponential integrators. I used both a one-stage and a seven-stage integrator. Both gave me the same hideous solution. It practically looks like noise.

I'm interested in using the HHT-$$\alpha$$ integrator, but I don't understand how to apply it to my system. Specifically, I was interested in using this paper. Could anyone explain it to me? How do I find the $$M$$ and $$Q$$ matrices? They have several distinct Jacobians, which one do I find if I take the Jacobian of the right-hand side of equation 4? How do I find the other Jacobian(s)?

When they perform the stability analysis on the simple equation $$u_t=\lambda u$$ I don't see how they even obtain the resultant equations. Considering this integrator was designed for rigid body mechanics I'm having difficulty translating the method to a more generic system of ODE.

And I know two of the papers are behind paywalls. I have the pdfs for both. Am I allowed to share them?

Does anyone know of some other integrators that might be useful?

Some properties that might be useful to know: The real and imaginary solutions constantly oscillate from positive to negative in time. It can be shown there is no steady-state solution for $$u$$. However, the term $$|u|^2$$ does reach a steady state.

MAJOR EDIT:

The full set of equations can be found in the paper "Effects of excitation inhomogeneity in semiconductor lasers pumped by an electron beam." As far as I can tell, this paper simply does not exist in PDF form on the internet. I had to scan a copy from my school's library. The equations to study are

$$\frac{\partial E}{\partial t}=\frac{\omega}{2}[\frac{i}{k^2}\frac{\partial ^2 E}{\partial x^2}+\frac{(1-i\chi_1)\epsilon ^"-\epsilon ^"_0}{\epsilon_r}E \tag{7}]$$ $$\frac{\partial \epsilon^"}{\partial t}=\frac{\sigma}{k}g_0g(x)-(\frac{1}{\tau_1}+2\sigma I)\epsilon^" \tag{8}$$ $$I=\chi_0 |E|^2 \tag{9}$$

making the change of variables to nondimensionalize the problem, I let $$x=\xi/k$$, $$t=s\tau_1$$, and $$E=u_0U$$ ($$u_0=1 V/m$$). This changes equations 7 and 8 to

$$\frac{\partial U}{\partial s}=\gamma(i\frac{\partial^2 U}{\partial \xi^2}+(\alpha \epsilon^"-\beta)U) \tag{10}$$ $$\frac{\partial \epsilon^"}{\partial s}=G_0g(\xi/k)-(1+\delta|U|^2)\epsilon^" \tag{11}$$

Where

$$\gamma=\frac{\tau_1\omega}{2}=1.07757*10^7$$

$$\alpha=(1-i\chi_1)/\epsilon_r=0.0775194 - 2.86484*10^{-6}*i$$

$$\beta=\epsilon^"_0/\epsilon_r=1.73648*10^{-7}$$

$$G_0=\tau_1\frac{\sigma}{k}*g_0=3.56747*10^6$$

$$\delta=2\tau_1\sigma\chi_0=0.0000111495$$

$$k=258195.$$

And the function g(x) is

$$g(x)=a*\exp[-b^2(x-\mu)^2]Sech(cx)$$

The coefficients are $$a=465.076$$, $$b=824.326$$, $$\mu=5.629393×10-3$$, $$c=-5633.573$$.

The domain goes from $$x=0$$ to $$x=L=6*10^{-3}$$.

The first thing to notice, is that if you use backwards euler on equation (11) you can find an exact expression (I'm dropping the double primes).

$$\epsilon^{n+1}_l=(\epsilon^n_l+\Delta sPG_0g(\xi_l/k))/(1+\Delta s(1+\delta |U^{n+1}_l|^2)) \tag{12}$$

Where I've added a new term P. This term is varied. At some values the solutions converge to zero, above a threshold, they (should) converge to a proper solution. For me, threshold seems to be about 2*10^-7.

If you apply the method I used in equation (6) you'll find it will converge with any initial condition you throw at it. For instance, setting the internal points to $$U(\xi,0)=\exp(i\xi)$$

• How are you calculating the Jacobian? A quick test in DifferentialEquations.jl seems to show that Rosenbrock23, Rodas5, and KenCarp4 do fine, but you just have to take care with the complex Jacobian (I just used the automatic analytical Jacobian computation). So I don't know your f but it seems to just work. I would make use of an integrator that's well-tested and adaptive, since the issue could be part implementation and non-adaptive methods generally fail on stiff equations. – Chris Rackauckas Jul 17 '19 at 10:24
• I'm using the standard method dF/du which gives J=iD+[df/du], where [] is a diagonal matrix with the derivative of f wrt u as its elements. I have been writing my methods in mathematica and I've been taking my matricies to be quite large (100). – Denis Jul 17 '19 at 11:55
• Your equation is very similar to the time-dependent Schrodinger equation, whose numerical solution has been studied extensively. I recommend using methods from that literature. – David Ketcheson Jul 17 '19 at 14:45
• Or maybe you are actually solving the nonlinear Schrodinger equation? That is also well-studied. You may get better answers if you reveal what $f$ actually is. – David Ketcheson Jul 17 '19 at 14:46
• Have you tried FDTD (Finite-Difference Time-Domain )? This method is heavily optimized for solving Schrodinger and Maxwell equations that usually involve imaginary variables or solutions. There are pretty good FDTD solvers out there such as this one: github.com/NanoComp/meep if you want to stick to Python. – Alone Programmer Nov 16 '19 at 18:52