# Improving calculation algorithm for coupled PDEs

I have the following two PDEs: $$\partial_zU=\nabla_r^2U+\varrho U$$ $$\partial_t\varrho=a\vert U\vert^4$$ with $a$ a constant and $$dt=dz\cdot\frac{n}{c}$$ with $n$ the refractive index of a material, and $c$ the speed of light.
My current approach is that I apply a crank-nicholson-algorithm on the first equation by neglecting the second part on the right side, using a matrix solver: $$\left(1-\frac{dz}{2.0}\partial_z\right)U_{n+1}=\left(1+\frac{dz}{2.0}\nabla_r^2\right)U_n$$ and then calculating $\varrho$ in between: $$U_{n+\frac{1}{2}}=\frac{U_{n+1}+U_{n}}{2}$$ $$\varrho_{n+\frac{1}{2}} = a\left\vert U_{n+\frac{1}{2}}\right\vert\cdot \underbrace{dz\cdot\frac{n}{c}}_{\equiv dt}$$ But is there an easier/more accurate way to do that, maybe even include it in the matrix equation above?

• Is U somehow time dependent? With your relation to z, you could rewrite it entirely with z-dependency only. $\partial_z(\frac{1}{U}(\partial_z U-\nabla^2_r U))=c/n a U^4$ and solve the nonlinear equation. This of course only works if U cannot be 0.
– Bort
Jul 13 '17 at 12:25
• Unfortunately U oscillates between positive and negative values, thus I expect it to be zero sometimes. Jul 13 '17 at 12:36

Given the dependency between $z$ and $t$, we can rewrite the first equation with a time derivative: $$\partial_tU=\frac{c}{n}\left(\nabla_r^2U+\rho U\right)$$. Now we define $g=\left(\begin{array}{c}U\\\rho\end{array}\right)$. This allows us to rewrite the system of equations as: $$\partial_t g=\left(\begin{array}{cc}\frac{c}{n}\nabla_r^2&0\\0&0\end{array}\right)g+\frac{c}{n}\left(\begin{array}{c}\frac{1}{2}\\0\end{array}\right)g^T\left(\begin{array}{cc}0&1\\1&0\\\end{array}\right)g+\left(\begin{array}{c}0\\1\end{array}\right)a\left(g^T\left(\begin{array}{cc}1&0\\0&0\end{array}\right)g\right)^2$$

Now you can choose your favorite time stepping algorithm for the non linear system of equations: $$\partial_tg=F(g,t,r)$$

E.g. Crank Nicolson would be a possible choice. Now there is no need to neglect something. Beware that you need to solve the non linearity for the implicit part of Crank Nicolson each step. For details see this question.

Edit

For the additional term $$\frac{c}{n} \beta|U|^2U$$. you need to add $$\frac{c}{n} \beta \left(\begin{array}{cc}1&0\\0&0\end{array}\right)g \left(g^T\left(\begin{array}{cc}1&0\\0&0\end{array}\right)g\right)$$. $U$ should not be in any expression, we want to solve for the field $g$.

For sure one can simplify the expressions further, but this should suffice as a starting point.

If you have more terms in $\partial_tU$, simply guess for $g$ (it's not that hard), evaluate it and see if you arrive at the original equations.

• Assumed I add $+\beta\vert U\vert^2U$ to the first equation, how would that part change? Am I correct if I simply add $+\frac{\beta n}{c}\begin{pmatrix}U\\0\end{pmatrix}\left(g^T\begin{pmatrix}1&0\\0&0\end{pmatrix}g\right)$ to the equation? Jul 13 '17 at 15:15
• Concerning the edit: After $\begin{pmatrix}1&0\\0&0\end{pmatrix}g=\begin{pmatrix}U\\0\end{pmatrix}$, there is no difference I can see? Thanks! Jul 13 '17 at 16:09
• Afaik the prefactor should be $\frac{c}{n}$, and not $\frac{n}{c}$. Can you confirm that? Jul 14 '17 at 14:38
• yes, it should be $\frac{c}{n}$. Let me fix the answer.
– Bort
Jul 14 '17 at 15:20

What you are doing is a kind of leap-frog scheme, but any time-stepping scheme will do. But since your problem is stiff and nonlinear, you may want to use a B-stable method, in order to avoid having to use small time steps for stability reasons.

• Would that approach also solve the problem that it is highly dependant on the size of z? In the current algorithm I get different results for different values of dz (which should not be the case). Jun 12 '17 at 14:58
• The values should converge if you choose $dx$ smaller and smaller. Jun 14 '17 at 14:35
• Unfortunately at least in my calculations it changes the result quite significantly when changing dz Jun 14 '17 at 14:55
• Can you set $n/c=1$ and see if things improve for small $dz=dt$? Jun 18 '17 at 12:43
• No, unfortunately not (except that the calculation time goes up when decreasing dz...) Jun 22 '17 at 7:45