- #1
TheCanadian
- 361
- 13
I am trying to solve a system of equations and have a question regarding the validity of my approach when implementing a fifth-order Cash-Karp Runge-Kutta (CKRK) embedded method with the method of lines. To give the questions some context, let me state the problem I am attempting to solve:
$$ \frac {\partial E}{\partial z} = - \frac {1}{c^2}\frac {\partial E}{\partial t} - \frac{1}{k} \frac {\partial ^2 E}{\partial z^2} - \frac{1}{kc^2} \frac {\partial^2 E}{\partial t^2} + iP \tag{1}
$$ $$
\\
\frac {\partial P}{\partial t} = iNE^* \tag{2}\\
$$ $$
\frac {\partial N}{\partial t} = iPE \tag{3}
$$
$$ E(z=0) = \frac{\partial E}{\partial z}(z=0) = E(t=0) = \frac{\partial E}{\partial t}(t=0) = 0,\\ P(t=0) = P_0e^{z/c}, N(t=0) = N_0e^{z/c} $$
where ## c = 3 \times 10^8, k = 1000, P_0 ## and ##N_0## are constants##, i=\sqrt{-1}####; 0 \leq t \leq 1000, 0 \leq z \leq 1000 ##
I am implementing CKRK on the above, and even though the first spatial derivative of ##E## depends on the second spatial derivative of ##E##, the numerical method appears to work when solving (1)-(3) when I use the scheme of approximating the time and spatial derivatives of ##E## on the right hand side of (1) by a backward difference approximation (I am using an accuracy of 5).
To switch (1) above to a system of first order spatial derivatives in ##z##, I could make the substitution:
$$ U = \frac {\partial E}{\partial z} $$
and solve the following equations for ##E## instead:
$$
U = \frac {\partial E}{\partial z} \tag{4}\\
$$ $$\frac {\partial U}{\partial z} = - \frac {k}{c^2}\frac {\partial E}{\partial t} - kU - \frac{1}{c^2}\frac {\partial^2 E}{\partial t^2} + kP \tag{5}
$$
But when testing these same initial/boundary conditions using the same numerical method on the coupled equations (2) - (5), the code takes too long to finish (the step sizes required become extremely small). I believe it is due to the fact that the coefficients on the right hand side of (5) are very large and cause stability issues. I have tried to rescale the values for ##z,t,P,N,## and ##E##, but doing so causes one of the other coupled equations to become unstable or has no effect (e.g. scaling ##z## does nothing to the value ## E = U\Delta z ## since both ##U## and ##\Delta z## would scale reciprocally and cancel any effect). It is due to similar reasons I am solving ##E## in the ##z##-direction as opposed to doing the substitution ##U = \frac {\partial E}{\partial t}## and solving it in ##t## which is the standard method of lines approach (when I tried this method, the ##\Delta t## given by CKRK becomes very small).
So ultimately, instead of using equations (2) - (5), I was wondering if applying CKRK to (1) - (3) is still a valid approach where I approximate the derivatives of ##E## on the right-side of (1) by backward finite differences? It seems very odd to apply CKRK to a first order spatial derivative that depends on an approximation of the second order spatial derivative, but is this wrong? (I would be using stored intermediate values of ##E## to ensure the backward finite difference approximations are also following the Runge-Kutta method.)
$$ \frac {\partial E}{\partial z} = - \frac {1}{c^2}\frac {\partial E}{\partial t} - \frac{1}{k} \frac {\partial ^2 E}{\partial z^2} - \frac{1}{kc^2} \frac {\partial^2 E}{\partial t^2} + iP \tag{1}
$$ $$
\\
\frac {\partial P}{\partial t} = iNE^* \tag{2}\\
$$ $$
\frac {\partial N}{\partial t} = iPE \tag{3}
$$
$$ E(z=0) = \frac{\partial E}{\partial z}(z=0) = E(t=0) = \frac{\partial E}{\partial t}(t=0) = 0,\\ P(t=0) = P_0e^{z/c}, N(t=0) = N_0e^{z/c} $$
where ## c = 3 \times 10^8, k = 1000, P_0 ## and ##N_0## are constants##, i=\sqrt{-1}####; 0 \leq t \leq 1000, 0 \leq z \leq 1000 ##
I am implementing CKRK on the above, and even though the first spatial derivative of ##E## depends on the second spatial derivative of ##E##, the numerical method appears to work when solving (1)-(3) when I use the scheme of approximating the time and spatial derivatives of ##E## on the right hand side of (1) by a backward difference approximation (I am using an accuracy of 5).
To switch (1) above to a system of first order spatial derivatives in ##z##, I could make the substitution:
$$ U = \frac {\partial E}{\partial z} $$
and solve the following equations for ##E## instead:
$$
U = \frac {\partial E}{\partial z} \tag{4}\\
$$ $$\frac {\partial U}{\partial z} = - \frac {k}{c^2}\frac {\partial E}{\partial t} - kU - \frac{1}{c^2}\frac {\partial^2 E}{\partial t^2} + kP \tag{5}
$$
But when testing these same initial/boundary conditions using the same numerical method on the coupled equations (2) - (5), the code takes too long to finish (the step sizes required become extremely small). I believe it is due to the fact that the coefficients on the right hand side of (5) are very large and cause stability issues. I have tried to rescale the values for ##z,t,P,N,## and ##E##, but doing so causes one of the other coupled equations to become unstable or has no effect (e.g. scaling ##z## does nothing to the value ## E = U\Delta z ## since both ##U## and ##\Delta z## would scale reciprocally and cancel any effect). It is due to similar reasons I am solving ##E## in the ##z##-direction as opposed to doing the substitution ##U = \frac {\partial E}{\partial t}## and solving it in ##t## which is the standard method of lines approach (when I tried this method, the ##\Delta t## given by CKRK becomes very small).
So ultimately, instead of using equations (2) - (5), I was wondering if applying CKRK to (1) - (3) is still a valid approach where I approximate the derivatives of ##E## on the right-side of (1) by backward finite differences? It seems very odd to apply CKRK to a first order spatial derivative that depends on an approximation of the second order spatial derivative, but is this wrong? (I would be using stored intermediate values of ##E## to ensure the backward finite difference approximations are also following the Runge-Kutta method.)