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