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?
2 Answers
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.
-
$\begingroup$ 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? $\endgroup$ Commented Jul 13, 2017 at 15:15
-
$\begingroup$ 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! $\endgroup$ Commented Jul 13, 2017 at 16:09
-
$\begingroup$ Afaik the prefactor should be $\frac{c}{n}$, and not $\frac{n}{c}$. Can you confirm that? $\endgroup$ Commented Jul 14, 2017 at 14:38
-
$\begingroup$ yes, it should be $\frac{c}{n}$. Let me fix the answer. $\endgroup$– BortCommented Jul 14, 2017 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.
-
$\begingroup$ 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). $\endgroup$ Commented Jun 12, 2017 at 14:58
-
$\begingroup$ The values should converge if you choose $dx$ smaller and smaller. $\endgroup$ Commented Jun 14, 2017 at 14:35
-
$\begingroup$ Unfortunately at least in my calculations it changes the result quite significantly when changing
dz$\endgroup$ Commented Jun 14, 2017 at 14:55 -
$\begingroup$ Can you set $n/c=1$ and see if things improve for small $dz=dt$? $\endgroup$ Commented Jun 18, 2017 at 12:43
-
$\begingroup$ No, unfortunately not (except that the calculation time goes up when decreasing
dz...) $\endgroup$ Commented Jun 22, 2017 at 7:45
Uoscillates between positive and negative values, thus I expect it to be zero sometimes. $\endgroup$