# 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.

• 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 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 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 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 at 14:46