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]);
B
? $\endgroup$ – Kirill Jul 31 '15 at 15:11