0

I am applying a 4th order Runge-Kutta code, using the method of lines, to solve the following:

\begin{equation} \frac {\partial y_1}{\partial t} = y_2 y_3 - C_1 y_1 \end{equation}

\begin{equation} \frac {\partial y_2}{\partial t} = y_3y_1 - C_2 y_2 \end{equation}

\begin{equation} \frac {\partial y_3}{\partial z} = y_4 \end{equation}

\begin{equation} \frac {\partial y_4}{\partial z} = \frac {\partial ^2 y_2}{\partial t^2} + \frac {\partial ^2 y_3}{\partial t^2} + y_1y_3 - y_2y_3 + y_3^2y_2 - y_4\end{equation}

(where $ y_n$ are the solutions I am looking for, $ C_n $ are real constants)

My code appears to work and solve the system when I just remove the $ y_4 $ term from equation (4), but otherwise, the solution goes to infinity. I have tried rescaling parameters/variables, decreasing the spatial intervals, and changing the boundary/initial conditions (except trivial ones), but these seem not to work. For testing purposes, I also changed equation (4) to:

$$ \frac {\partial y_4}{\partial z} = \frac {\partial ^2 y_2}{\partial t^2} + \frac {\partial ^2 y_3}{\partial t^2} + y_1y_3 - y_2y_3 + y_3^2y_2 - \text {SCALE}*y_4 $$

where SCALE = $10^{-5}$. With this change, the output did not go to infinity, although it had an unexpected output (namely there were frequent large oscillations). I could not change SCALE to values larger than $10^{-5}$, otherwise the output would go to infinity.

RK4 seems to work well for almost everything except this one term in the 4th equation, and I was just wondering if you had any ideas? I was considering possibly applying an implicit method on only the 3rd and 4th equations, although due to the nonlinearity, I haven't been able to isolate for the $y_3$ when using higher order BDF methods.

Mathews24
  • 578
  • 2
  • 6
  • 13
  • 2
    Runge-Kutta methods solve systems of first-order ODEs. That is not what you have here. You need to rewrite it as a first-order system. – David Ketcheson Dec 10 '16 at 18:25

1 Answers1

-2

This is definitely a stiffness issue. Instead of using RK4, what about using a reasonably well-optimized adaptive timestepping method? In many cases that will work well enough. If not, the places to go are Rosenbrock methods are BDF methods. You don't need to isolate anything to use BDF methods: you know the second time derivatives (since you know the first derivatives), so you can write down the vector function for how everything changes and just plop that into an ODE solver (with all 4 equations). The ODE solver will solve the root finding problem for you.

Chris Rackauckas
  • 12,320
  • 1
  • 42
  • 68
  • Thank you for the response. I tried an adaptive tilmestep using the Fehlberg method but it was going too slowly as there were small oscillations present and the time steps were on the order of nanoseconds while my domain is days. I will look into implicit methods, but I don't have an explicit expression for $\frac {\partial ^2 y_2}{\partial t^2}$ (nor the first derivative in time), so will it work? – Mathews24 Oct 09 '16 at 16:13
  • The same expression that you're plugging into the ODE solver for the RK method should just plug into the BDF method. Also, I don't know what ODE solver you're using but you shouldn't use a Fehlberg method, those are really awful in comparison to modern methods like Dormand-Prince and are never suggested. – Chris Rackauckas Oct 09 '16 at 16:17
  • Hmmm, interesting. I looked back and it was actually Cash-Karp that I used, but isn't the difference just 4th order vs 5th order? Also, I was under the impression that although lower order methods were less accurate, they had a greater radius of convergence? Maybe I am simply lacking the understanding, but my equation for $y_4$ is a partial differential equation with respect to $z$ (and my equation for $y_3$ is a PDE in $z$), where as BDF methods would require them all to be differential equations of one variable (e.g. $t$), no? – Mathews24 Oct 09 '16 at 20:24
  • That can be true. Sometimes a higher order method can have the same or a larger radius of convergence. Cash-Karp and RKF are both order 4/5, but they have much different properties. See this post. Any ODE solver you use will propogate it forward in terms of the derivatives you specify it, i.e. if you natively throw this into an ODE solver it will go forward in time for $y_1$ and $y_2$, and in $z$ for $y_3$ and $y_4$. – Chris Rackauckas Oct 09 '16 at 20:30
  • But what I am saying is, if you made a function f to put into one ODE solver, there's no reason why that same f cannot be used with another method. – Chris Rackauckas Oct 09 '16 at 20:31
  • Thank you for the comments. I've been writing my own code to learn more about these so haven't yet used an established ODE package yet, but I'll certainly check them out. I'll also continue looking into the Cash-Karp method I tried implementing earlier. I guess my other concern was wondering if the stiff ODE packages would provide an appropriate adaptive step if there were multiple derivative types given, but I'll have to check that out. If you have any ideas on packages for Python and implementation tips, I'd be glad to read them; otherwise, I'll get back to you (hopefully soon) with results. – Mathews24 Oct 09 '16 at 20:43