Finite Difference Beam Propagation Method problem

Question

I am trying to implement the finite difference beam propagation method to study the propagation of a TE light signal through a waveguide. However, my solutions are exponentially growing, and display no propagation. For more information, see "Introduction to Optical Waveguide Analysis" by Kawano and Kitoh, section 4.3.

This method involves iteratively solving a complex-valued matrix equation of the form $Ax = b$, with $A$ a square matrix, $b$ a known vector and $x$ a vector to be solved for. I am using the LAPACK routine zgtsv via C++ to solve this equation, but I am getting numerical instabilities. In particular, the vectors are exponentially growing.

More explicitly, given the complex function $\phi(x,0)$ we wish to solve for $\phi(x,z)$ subject to the partial differential equation

$$ 2i\beta\frac{\partial \phi}{\partial z} = \frac{\partial^2 \phi}{\partial x^2}+(k_0^2 n(x)^2 - \beta^2)\phi$$

where $-3<x<3$ and $k_0=4$. $n(x) = 3.5$ for $-1<x<1$ and $n(x) = 1.5$ otherwise. Suppose further that $0<z<40$. $\beta = 14$ is the propagation constant of the waveguide.

We discretize the problem into $N$ steps in the $x$ direction and $M$ steps in the $z$ direction and replace $\phi(x,z)$ by $\phi(p,q)$ where $p$ and $q$ are integers. Let $ dx = \frac{6}{N-1} $ and $ dz = \frac{40}{M-1} $ Using $(p,q+1/2)$ as a center difference for Taylor expansions, we get the equation

$$-\frac{1}{dx^2} \phi(p-1,q+1) + \left\{\frac{2}{dx^2} + \frac{4i\beta}{dz} - (k_0^2 n(x)^2 - \beta^2)\right\}\phi(p,q+1)-\frac{1}{dx^2} \phi(p+1,q+1) = \frac{1}{dx^2} \phi(p-1,q) + \left\{-\frac{2}{dx^2} + \frac{4i\beta}{dz} + (k_0^2 n(x)^2 - \beta^2)\right\}\phi(p,q)+\frac{1}{dx^2} \phi(p+1,q)$$ which much be satisfied for all $p$. This allows us to step forward from $\phi(p,q)$ to $\phi(p,q+1)$. This equation may be rephrased as a matrix equation $A \phi(q) = b$, where $A$ is a symmetric tridiagonal matrix.

However, this method does not appear to be working for me. In fact, $\phi(p,q)$ is growing exponentially with $q$.

Are there any known problems with using zgtsv iteratively like this? I have very little idea of what is going wrong.

EDIT: Here is a condensed version of the code I am using. I am using Dislin to plot Solution.

double PI = 3.141592;
std::complex<double> imaginaryI(0.0,1.0);

int N = 300;
int M = 300;

std::complex<double> *ComplexSolution= new std::complex<double> [M*N];

for (int i = 0; i < N; i++){
    std::cout << i << std::endl;
    if( i == N /2)
        ComplexSolution[i] = 1.0;
    else ComplexSolution[i] = 0.0;}

double *Solution= new double [M*N];
double waveK = 4.0;
double x1 = -3.0, x2 = -1.0, x3 = 1.0, x4 = 3.0;
double indexN1 = 1.5, indexN2 = 3.5;
double Length = 40;
double dx = (x4-x1)/(N-1), dz = Length/(M-1);
double BETA = 14.0;

double alpha1 = 1.0/(dx*dx), alpha2 = -2.0/(dx*dx);

std::complex<double> *BPM_sub_diagonal = new std::complex<double> [N-1];
std::complex<double> *BPM_diagonal = new std::complex<double> [N];
std::complex<double> *workBPM_sub_diagonal = new std::complex<double> [N-1];
std::complex<double> *workBPM_diagonal = new std::complex<double> [N];
std::complex<double> *B = new std::complex<double> [N];

int NRHS = 1;
int INFO;

for (int i = 0; i < N-1; i++)
    BPM_sub_diagonal[i] = -alpha1;

for (int i = 0; i < N; i++){
    if(i*dx+x1>x2 && i*dx+x1<x3)
        BPM_diagonal[i] = -alpha2 + (4.0*imaginaryI*BETA/dz) - waveK*waveK*indexN2*indexN2 + BETA*BETA;
    else
        BPM_diagonal[i] = -alpha2 + (4.0*imaginaryI*BETA/dz) - waveK*waveK*indexN1*indexN1 + BETA*BETA;}

for(int j = 1; j < M; j++){
    for(int i = 0; i < N; i++)
        workBPM_diagonal[i] = BPM_diagonal[i];

    for(int i = 0; i < N-1; i++)
        workBPM_sub_diagonal[i] = BPM_sub_diagonal[i];

    for(int i = 0; i < N; i++){
        if(i == 0)
            B[i] = (alpha2 + (4.0*imaginaryI*BETA/dz) + waveK*waveK*indexN1*indexN1 - BETA*BETA)*ComplexSolution[(j-1)*N+i]
                   + alpha1*ComplexSolution[(j-1)*N+i+1];

        else if(i == N-1)
            B[i] = (alpha2 + (4.0*imaginaryI*BETA/dz) + waveK*waveK*indexN1*indexN1 - BETA*BETA)*ComplexSolution[(j-1)*N+i]
                   + alpha1*ComplexSolution[(j-1)*N+i-1];

        else if(i*dx+x1>x2 && i*dx+x1<x3)
            B[i] = (alpha2 + (4.0*imaginaryI*BETA/dz) + waveK*waveK*indexN2*indexN2 - BETA*BETA)*ComplexSolution[(j-1)*N+i]
                   + alpha1*ComplexSolution[(j-1)*N+i+1]
                   + alpha1*ComplexSolution[(j-1)*N+i-1];
        else 
            B[i] = (alpha2 + (4.0*imaginaryI*BETA/dz) + waveK*waveK*indexN1*indexN1 - BETA*BETA)*ComplexSolution[(j-1)*N+i]
                   + alpha1*ComplexSolution[(j-1)*N+i+1]
                   + alpha1*ComplexSolution[(j-1)*N+i-1];
    }

    zgtsv_(&N, &NRHS, workBPM_sub_diagonal, workBPM_diagonal, workBPM_sub_diagonal, B, &N, &INFO);

    for(int i = 0; i < N; i++)
        ComplexSolution[j*N+i] = B[i];

}

for(int i = 0; i < M*N; i++)
    Solution[i] = std::abs(ComplexSolution[i]);

Answer 1

When you call Lapack's zgtsv, it doesn't just solve a tridiagonal system $Ax=b$. What it does first is perform an LU factorization (zgttrf) $A = LU$, where $L,U$ are lower- and upper-tridiagonal matrices, and only then proceeds to solve $LUx=b$.

When you give it the lower, main, and upper diagonals of the matrix $A$, those diagonals are overwritten by the routines that zgtsv invokes to contain, on output, the lower diagonal of $L$, main diagonal of $U$ and upper diagonal of $U$, respectively. This is a standard idea in linear algebra libraries like Blas and Lapack.

This only needs to be done once per matrix: the LU factorization can be reused for future solves of $Ax=b$ (zgttrs).

Your code does the following thing wrong: it uses the same storage for the lower diagonal of $A$ and for its upper diagonal, meaning that Lapack will necessarily clobber its own output when writing the matrices $L$ and $U$ into the storage of $A$. The solution is to introduce another working vector

workBPM_super_diagonal[i] = BPM_sub_diagonal[i];

and pass that to zgtsv instead of passing workBPM_sub_diagonal to it twice.

Here's the iteration $j=10$ when I did this:

I am not too sure what results you are expecting there, but it doesn't blow up (for $N,M=300$). Since Crank-Nicolson wouldn't do a good job of damping oscillations due to discontinuous coefficients or delta-function initial functions, I'm not sure whether all of those oscillations are expected from the equation or if some of them are spurious.

Other things:

• Reuse the LU factorization between time steps, as your matrix does not change.
• Do not use raw pointers: it is not necessary given that there are standard C++ types such as vector and valarray. Learn about standard C++ features that C++ expects you to use.
• See The Definitive C++ Book Guide and List: https://stackoverflow.com/q/388242/491171
• As a principle, avoid duplication: the r.h.s. of your equation is specified in six different places. If there is an error, or if you want to solve an equation with a different r.h.s., that is five places too many that need to be correct at the same time.
• Decouple your code into modular functions that each solve a small, carefully-defined task. There is no reason why the definition of your equation has to be included right in the FD solver. Unnecessary coupling means that different parts of your program cannot be tested in isolation. C++ includes many tools to design modular testable programs.
• Use M_PI.
• Consider asking on https://codereview.stackexchange.com/.