Next: An example solution of Up: Poisson's equation Previous: 1-d problem with mixed

## An example 1-d Poisson solving routine

Listed below is an example 1-d Poisson solving routine which utilizes the previously listed tridiagonal matrix solver and the Blitz++ library (see Sect. 2.20).
```// Poisson1D.cpp

// Function to solve Poisson's equation in 1-d:

//  d^2 u / dx^2 = v  for  xl <= x <= xh

//  alpha_l u + beta_l du/dx = gamma_l  at x=xl

//  alpha_h u + beta_h du/dx = gamma_h  at x=xh

// Arrays u and v assumed to be of extent N+2.

// Now, ith elements of arrays correspond to

//  x_i = xl + i * dx    i=0,N+1

// Here, dx = (xh - xl) / (N+1) is grid spacing.

#include <blitz/array.h>

using namespace blitz;

void Tridiagonal (Array<double,1> a, Array<double,1> b, Array<double,1> c,
Array<double,1> w, Array<double,1>& u);

void Poisson1D (Array<double,1>& u, Array<double,1> v,
double alpha_l, double beta_l, double gamma_l,
double alpha_h, double beta_h, double gamma_h,
double dx)
{
// Find N. Declare local arrays.
int N = u.extent(0) - 2;
Array<double,1> a(N+2), b(N+2), c(N+2), w(N+2);

// Initialize tridiagonal matrix
for (int i = 2; i <= N; i++) a(i) = 1.;
for (int i = 1; i <= N; i++) b(i) = -2.;
b(1) -= beta_l / (alpha_l * dx - beta_l);
b(N) += beta_h / (alpha_h * dx + beta_h);
for (int i = 1; i <= N-1; i++) c(i) = 1.;

// Initialize right-hand side vector
for (int i = 1; i <= N; i++)
w(i) = v(i) * dx * dx;
w(1) -= gamma_l * dx / (alpha_l * dx - beta_l);
w(N) -= gamma_h * dx / (alpha_h * dx + beta_h);

// Invert tridiagonal matrix equation
Tridiagonal (a, b, c, w, u);

// Calculate i=0 and i=N+1 values
u(0) = (gamma_l * dx - beta_l * u(1)) /
(alpha_l * dx - beta_l);
u(N+1) = (gamma_h * dx + beta_h * u(N)) /
(alpha_h * dx + beta_h);
}
```

Richard Fitzpatrick 2006-03-29