Python Runge Kutta for nonlinear system of equation

[TheCanadian]

I am applying a 4th order Runge-Kutta code 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.

 

[DrClaude]

You have dependences in z and t. That means PDEs, not ODEs. How to you solve that with RK?

 

[TheCanadian]

You have dependences in z and t. That means PDEs, not ODEs. How to you solve that with RK?

I was using the method of lines, advancing $y_1$ and $y_2$ first, then computing the values for $y_3$ and $y_4$.

 

[DrClaude]

I was using the method of lines, advancing $y_1$ and $y_2$ first, then computing the values for $y_3$ and $y_4$.

Ok. Have you looked up the criterion of stability for that approach?

 

[TheCanadian]

Ok. Have you looked up the criterion of stability for that approach?

I had previously looked into stability conditions (e.g. Von Neumann) although it appeared there were no standard criterion for dealing with nonlinear systems that presented stiff behaviour.

 

[ChrisVer]

what kind of system does these equations describe? Or they are just made up?
The fact that introducing a small scale didn't make your solutions explode, at least indicates that something is stinks with that term... did you try higher values of SCALE to see how things change?

 

[TheCanadian]

what kind of system does these equations describe? Or they are just made up?
The fact that introducing a small scale didn't make your solutions explode, at least indicates that something is stinks with that term... did you try higher values of SCALE to see how things change?

There are a lot of real and imaginary constants I omitted, but the general form of the equations above models an actual physical system of coherently radiating molecules (i.e. it is not made up). Yes, when I increase SCALE, the values go to infinity. When I decrease scale, the large oscillations are gone and the output resembles what I expect for a particular case of initial/boundary conditions, but there are quite a few differences that I'm thinking are largely due to this one term being absent when SCALE approaches 0.