A tutorial on C++ based GSL code for solving partial differential equations (PDE) with finite differencing

Presented below is a C++ based code example of using GSL for solving partial differential equations with finite differencing.

The example problem is the free vibration case of an Euler Bernoulli beam.

Field equation:

EI (d4w/dx4) + (rho*A) (d2w/dt2) = 0

Boundary conditions:

Clamped:

w = 0, dw/dx = 0

Hinged:

w = 0, d2w/dx2 = 0

Free:

d3w/dx3 = 0, d2w/dx2 = 0

Notations are the usual ones.

The code begins

#include < iostream >
#include < fstream >
#include < math.h >

//REF.: http://www.gnu.org/software/gsl/manual/html_node/Linear-Algebra-Examples.html#Linear-Algebra-Examples
#include < stdio.h >
#include < gsl/gsl_linalg.h >
//REF.: http://www.gnu.org/software/gsl/manual/html_node/Linear-Algebra-Examples.html#Linear-Algebra-Examples

//REF.: https://www.gnu.org/software/gsl/manual/html_node/BLAS-Examples.html#BLAS-Examples
#include < gsl/gsl_eigen.h >
#include < gsl/gsl_blas.h >
//REF.: https://www.gnu.org/software/gsl/manual/html_node/BLAS-Examples.html#BLAS-Examples

using namespace std;

int main (){

 cout<<"Hello, i am eulerbeam1... \n";

 fstream ofl1;
 ofl1.open("eulerbeam1.res", ios::out | ios::ate);

 //REF.: http://www.techscience.com/doi/10.3970/sl.2011.006.115.pdf
 double E = 2.1e11, rho = 7.85e3, l = 1, b = 0.01, d = 0.005;
 double I = b*d*d*d/12, A = b*d;
 
 int Nnodes = 60;

 double x0 = 0, xn = 1;
 double deltax = (xn - x0)/(Nnodes - 1);
 double onb2deltax = 1./(2*deltax);

 gsl_vector * vector1 = gsl_vector_alloc(Nnodes);

 gsl_matrix * Cmatx_nt = gsl_matrix_calloc (Nnodes, Nnodes);//Coefficients for first derivative
 gsl_matrix * Cmatx_nt2 = gsl_matrix_calloc (Nnodes, Nnodes);//Coefficients for second derivative
 gsl_matrix * Cmatx_nt3 = gsl_matrix_calloc (Nnodes, Nnodes);//Coefficients for third derivative
 gsl_matrix * Cmatx_nt4 = gsl_matrix_calloc (Nnodes, Nnodes);//Coefficients for fourth derivative
 gsl_matrix * Cmatx_Id = gsl_matrix_calloc (Nnodes, Nnodes);//the Identitiy matrix
 gsl_matrix * Cmatx_K = gsl_matrix_calloc (Nnodes, Nnodes);//the Stiffness matrix
 gsl_matrix * Cmatx_M = gsl_matrix_calloc (Nnodes, Nnodes);//the Mass matrix
 gsl_matrix * Cmatx_invK_M = gsl_matrix_calloc (Nnodes, Nnodes);//for storingthe inv(K)*M

 //for 0th node
 gsl_matrix_set(Cmatx_nt, 0, 0, -1./deltax);
 gsl_matrix_set(Cmatx_nt, 0, 1, 1./deltax);
 //for nth node
 gsl_matrix_set(Cmatx_nt, Nnodes - 1, Nnodes - 2, -1./deltax);
 gsl_matrix_set(Cmatx_nt, Nnodes - 1, Nnodes - 1, 1./deltax);
 for(int i = 1; i < Nnodes - 1; i++){
  gsl_matrix_set(Cmatx_nt, i, i - 1, -onb2deltax);
  gsl_matrix_set(Cmatx_nt, i, i + 1, onb2deltax);
 }

 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
     1.0, Cmatx_nt, Cmatx_nt,
     0.0, Cmatx_nt2);
 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
     1.0, Cmatx_nt, Cmatx_nt2,
     0.0, Cmatx_nt3);
 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
     1.0, Cmatx_nt2, Cmatx_nt2,
     0.0, Cmatx_nt4);

 gsl_matrix_set_identity(Cmatx_Id);


 gsl_matrix_memcpy(Cmatx_K, Cmatx_nt4);

 gsl_matrix_get_row(vector1, Cmatx_Id, 0);//for (deflection, w = 0) at x = 0
 //gsl_matrix_get_row(vector1, Cmatx_nt3, 0);//for (shear force, d3w/dx3 = 0) at x = 0
 gsl_matrix_set_row(Cmatx_K, 0, vector1);

 gsl_matrix_get_row(vector1, Cmatx_Id, Nnodes - 1);//for (deflection, w = 0) at x = 1
 //gsl_matrix_get_row(vector1, Cmatx_nt3, Nnodes - 1);//for (shear force, d3w/dx3 = 0) at x = 1
 gsl_matrix_set_row(Cmatx_K, Nnodes - 1, vector1);

 //gsl_matrix_get_row(vector1, Cmatx_nt, 0);//for (slope, dw/dx = 0) at x = 0
 gsl_matrix_get_row(vector1, Cmatx_nt2, 0);//for (moment, d2w/dx2 = 0) at x = 0
 gsl_matrix_set_row(Cmatx_K, 1, vector1);

 //gsl_matrix_get_row(vector1, Cmatx_nt, Nnodes - 1);//for (slope, dw/dx = 0) at x = 1
 gsl_matrix_get_row(vector1, Cmatx_nt2, Nnodes - 1);//for (moment, d2w/dx2 = 0) at x = 1
 gsl_matrix_set_row(Cmatx_K, Nnodes - 2, vector1);

 gsl_matrix_set_identity(Cmatx_M);
 gsl_matrix_set(Cmatx_M, 0, 0, 0);
 gsl_matrix_set(Cmatx_M, 1, 1, 0);
 gsl_matrix_set(Cmatx_M, Nnodes - 2, Nnodes - 2, 0);
 gsl_matrix_set(Cmatx_M, Nnodes - 1, Nnodes - 1, 0);

//Computing inv(K)*M through LU decomposition
 int signum; 
 gsl_permutation * p_lu = gsl_permutation_alloc (Nnodes);
 gsl_linalg_LU_decomp(Cmatx_K, p_lu, &signum);
 for(int i=0; i < Nnodes; i++){
  gsl_matrix_get_col (vector1, Cmatx_M, i); 
  gsl_linalg_LU_svx(Cmatx_K, p_lu, vector1);
  gsl_matrix_set_col(Cmatx_invK_M, i, vector1);
 }


//Computing eigenvalues and eigenvectors
 gsl_vector_complex *eval = gsl_vector_complex_alloc (Nnodes);
 gsl_matrix_complex *evec = gsl_matrix_complex_alloc (Nnodes, Nnodes);
 gsl_eigen_nonsymmv_workspace * myeigen_ws =
  gsl_eigen_nonsymmv_alloc (Nnodes);
 gsl_eigen_nonsymmv (Cmatx_invK_M, eval, evec, myeigen_ws);
 gsl_eigen_nonsymmv_free (myeigen_ws);
 gsl_eigen_nonsymmv_sort (eval, evec,
   GSL_EIGEN_SORT_ABS_DESC);
  

//saving the frequencies(Hz) in ascending order
 int i, j;
 double dumi= (1./(2*M_PI*l*l))*sqrt(E*I/(rho*A));
 for (i = 0; i < Nnodes; i++)
 {
  gsl_complex eval_i
   = gsl_vector_complex_get (eval, i);
  gsl_vector_complex_view evec_i
   = gsl_matrix_complex_column (evec, i);
  ofl1 << sqrt(1./(GSL_REAL(eval_i)))*dumi << " i" << GSL_IMAG(eval_i) << endl;
 }



 ofl1.close();
 cout << "Hello, i was eulerbeam1... \n";
 return 0;
}

The code ends.
Compilation command on a Linux console/command line:

g++ -lgsl eulerbeam1.cc -o eulerbeam1

Execution:
./eulerbeam1

Results give the following fundamental frequencies (Ref.: http://www.techscience.com/doi/10.3970/sl.2011.006.115.pdf):

Clamped-Clamped: 26.5281Hz
Clamped-Free: 4.44362 Hz
Simply-Supported: 11.7253 Hz