# Finite Difference Beam Propagation Method problem

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]);

• Can you be more specific about what exactly you are solving, what the iteration is and how it's related to the linear system? I don't have that book, so I don't know what you're doing exactly, I suspect this would be true for other people as well. It's confusing as is, because zgtsv solves a tridiagonal system directly by GE w/pivoting, so it is unlikely to be the main cause of problems. What is $A$? What iteration? "Numerical instabilities" is also a bit too vague, can you expand? What vectors are growing, and how exactly? – Kirill Jul 30 '15 at 11:45
• Does the blow-up also happen when you increase $M$ significantly? Is it slower? Also, you probably already did this (several times), but you should check your implementation of $A$ for incorrect signs (watch out for wrong parentheses, too). Getting a single sign wrong could push the eigenvalues of $A$ into the instable region, and hence you get blow-up. – Christian Clason Jul 30 '15 at 14:49
• Is it $2i$ or $2i\beta$? I think the error could be in incorrect implementation of $A$ (can't be sure without seeing code), but it can't be the zgtsv itself causing problems, as all it does is solve the linear equations. The discretization (Crank-Nicolson) looks fine. Can you post the code? – Kirill Jul 30 '15 at 15:07
• @Kirill: Good catch! The Crank-Nicolson discretization looks a bit obfuscated here. – Geoff Oxberry Jul 31 '15 at 10:31
• Isn't the sign of $k_0^2n(x)^2-\beta^2$ wrong in the r.h.s. vector B? – Kirill Jul 31 '15 at 15:11

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.