# Crank–Nicolson method for nonlinear differential equation

I want to solve the following differential equation from a paper with the boundary condition:

The paper used the Crank–Nicolson method for solving it. I think I understand the method after googling it, but most websites discussing it use the heat equation as an example. Here we can replace the usual t variable with Xi, and the usual x as Rho. Now unlike the heat equation, this equation has an additional first derivative of Xi term (the first term gives both first and second derivative of Xi). Does it matter what finite difference scheme I use to approximate this term?

Another thing is that the equation has nonlinear terms involving |A|^2 and |A|^4. How can I incorporate this to the matrix equation I am going to solve? (PS: A is complex)

I am pretty new to numerical methods so forgive me if I have neglected anything significant.

• There are similar questions with Crank Nicolson and nonlinear PDEs, e.g. this, where I have explained the concept and how it is applied. What you need first is a solver for the linear problem and then use the Newton-Kantorovich scheme on every "time step" with the Frechet Derivative of your functional. – Bort Feb 12 '18 at 14:47

You need at least one initial condition for $\xi$ and two BC for $\rho$. Where are they? Your equation looks like the heat equation in cylindrical coordinates assuming angular and plane symmetry, with nonlinear heat sources: $$2ik\partial_tA-\triangle A=N(A)$$.

I would discretise the spatial derivative with finite differences (which is essentialy the same as finite elements in this case) and discretise the dynamic system with what you've proposed. If you discretise in space you will obtain: $$\partial_t\vec{A}-L(\rho) \vec{A}=\vec{N}(\vec{A})$$ Where $L(\rho)$ is the matrix operator for the Laplacian.

This operator looks like $$L(\rho)=L+RG$$ Where $L$ is the famous Laplacian matrix operator in cartesian coordinates and $G$ the gradient operator. Without any BC these operators look like: $$L(i,i)=-2/\Delta \rho^2, L(i-1,i)=L(i,i+1)=1/\Delta \rho^2$$ $$G(i,i+1)=-G(i,i-1)=1/\Delta\rho$$ $$R(i,i)=1/\rho_i$$

The resulting scheme will be sth. like this: $$[2ki\mathbb{I}-\Delta \xi/2L(\rho)]\vec{A}^{n+1}+\Delta\xi/2\vec{N}(\vec{A}^{n+1})= 2ki\vec{A}^n+\Delta\xi/2 L(\rho)\vec{A}^{n} +\Delta\xi/2\vec{N}(\vec{A}^{n})$$

Has been obtained from the general Crank-Nicolson scheme: $$\frac{W^{n+1}-W^n}{\Delta t}=\frac{1}{2}[F(W^{n+1})+F(W^n)]+\mathcal{O}(\Delta t^2)$$ for the differential equation: $$\frac{dW}{dt}=F(t,W)$$

For each $n$ you will solve a nonlinear system of equations...

• Using finite differences in cylindrical coordinate systems is difficult, because you need to do something with the node at $\rho=0$ where the equation is singular. It's simpler to do this with finite elements. – Wolfgang Bangerth Feb 6 '18 at 18:33
• Also, you really only need one boundary condition at $\rho=R$. No boundary condition is necessary at $\rho=0$ -- it is implied by the symmetry of the solution. – Wolfgang Bangerth Feb 6 '18 at 18:34
• The simplest method (in cartesian grid) is FD whatever the equation, in my opinion. Difficult would be if used for solving problems with curved computational domain. The singularity is avoidable if Neuman BC is provided, and it does not matter if Dirichlet, so there is not any problem with that issue. You are assuming that is a solid cylinder, I am considerig a more general case with a hole. If you compute the matrix operators you will find that actually two conditions must be given for the uniqueness of the solution. Name them as you wish. – HBR Feb 7 '18 at 11:04
• If it's a cylindrical shell, you are correct. If the cylinder is solid, mathematically you can't pose a Neumann-type boundary condition. – Wolfgang Bangerth Feb 7 '18 at 14:43
• So the second equation A_entrance is basically the initial condition for ξ at ξ=0 (not to be confused with z, which is just a parameter that I input). The equation describes the propagation of electric-field in nonlinear medium, so A tends to zero as rho tends to infinity as the A has a Gaussian profile initially (and shouldn't deviate from Gaussian shape too much after the whole propagation). Are there anymore BC/IC I need? – Physicist Feb 9 '18 at 9:13