Implementing odespy for system of PDEs

Question

After trying to use RK4 to solve the below system of equations, it appears the output had large and fast oscillations even with an adaptive time step I incorporated using the Cash-Karp method. I am now looking to existing solvers with more sophisticated tools (e.g. BDF methods) to solve PDEs, and came across odespy which I hope will work.

Here is the problem:

Tutorials for odespy are fairly detailed. Although helpful, an example of how to implement a system of PDEs with both time and spatial derivatives is missing, and I am wondering if odespy can actually solve such a problem. I am mainly curious since using the method of lines reduces this case to solving a system of a system of ODEs. Here is my code thus far:

import odespy
import matplotlib
matplotlib.use('pdf')
import os 
import matplotlib.pyplot as plt 
from pylab import *
import numpy as np 

c = 3.*(10.**8.)  
N = 500
zmax = 1000.
z = np.linspace(0.,zmax,N+1)
f_kwargs = dict(L=zmax, z=z)
y_1 = np.zeros((N+1),dtype=np.complex_) #should be [Nz x Nt] array
y_2 = np.zeros((N+1),dtype=np.complex_)
y_3 = np.zeros((N+1),dtype=np.complex_)
y_4 = np.zeros((N+1),dtype=np.complex_)
u = y_1,y_2,y_3,y_4

U_0 = np.zeros((4,N+1),dtype=np.complex_) 
U_0[0][:] = np.cos(0.01)*(np.e**(z[:]/c))
U_0[1][:] = np.sin(0.01)*(np.e**(z[:]/c))
U_0[2][:] = 0. 
U_0[3][:] = 0. 

def rhs(u, t, L=None, beta=None, z=None): 
    N = len((u[0])[:]) - 1 #the PDE example specifies end condition, uses -1
    dz = z[1] - z[0] 
    for i in range(0, N):
        if i == 0:
            OUT = [np.cos(0.01)*np.e**(-t),
                   np.sin(0.01)*np.e**(-t),
                    0.,
                    0.] 

        else:
            OUT = [(y_2[i]*np.conj(y_3[i]) - y_1[i]),
                   (np.conj(y_3[i]*y_1[i]) - y_2[i]),
                    y_4[i],
                (y_4[i] - ((y_3[i])**2.)*(1./dz)*np.conj(y_2[i+1] - y_2[i]))]
    return OUT

solver = odespy.ForwardEuler(rhs, f_kwargs=f_kwargs,complex_valued=True)
solver.set_initial_condition(U_0)#solver.set_initial_condition(U_0)

dz = z[1] - z[0]
time_points = np.linspace(0.,600.,600) #t' = t/TR
dt = time_points[1] - time_points[0]

u, t = solver.solve(time_points)

This code does not work, and my main concerns are that I am:

1) Unable to properly formulate the equations in such a way odespy could handle a time derivative in my PDE (the 4th equation), and that the 3rd and 4th equations were derivatives in z as opposed to t (which is the derivative of the first two equations).

2) Unable to provide input to odespy such that it could take a system of PDEs, and then apply the method of lines approach outlined here. For example, one error already visible is my initial/boundary condition implementation. I used a 4 x (N+1) array for my output and initial/boundary conditions, but the solver is not expecting such an input from what I've gathered.

I'm still going through the documentation, but was wondering if anyone had any input for the above? Any suggestions at all on how to solve the above problem and comments on using odespy/numerical methods for a system of PDEs would be greatly appreciated.

Answer 1

Odespy needs from you a 1st order system of ODEs, namely something of the form \begin{align} \frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y},t), \tag{1} \end{align}

where, for your particular problem, $y_i$ are single-variate complex-valued functions that depend solely on time (that means you need to arrive at the semi-discrete form first by discretizing the spatial dimension, this is what is meant by the method of lines), and $i=1 \ldots3 $. Note that $y_4$ is just a placeholder variable defined (for the original space-time problem) by $y_4 := \partial y_3 / \partial z$.

So if we get rid of $y_4$ in your original formulation, we get the space-time system: \begin{eqnarray} \frac{\partial y_1}{\partial t} &=& \mbox{(as before)} \\ \frac{\partial y_2}{\partial t} &=& \mbox{(as before)} \\ \frac{\partial y_3}{\partial t} &=& \frac{\partial^2 y_3}{\partial z^2} - \frac{\partial y_3}{\partial z} + y_3^2 \frac{\partial y_2^*}{\partial z}. \end{eqnarray}

So what's left to do is:

1. discretize the RHS in space ($z$-terms) using suitable finite-difference formulas for the second and first partial derivatives in $z$, to arrive at a semi-discrete system like (1) above, and
2. pass the semi-discrete system to Odespy (make sure you get the boundary conditions right).

Make sure to choose a time integration method that can handle the time scales and other properties of your system (for example, forward Euler does not work that well for stiff systems).

Later edit: had a closer look at the system and the $y_i$ look to be real-valued functions of real variables (initial/boundary conditions and right-hand side are all real). So if you know your solutions have nonzero imaginary parts, you should double-check your PDE system.