Solves the harmonic equation. More...

#include <Utility.h>
#include <EVP_bundle.h>
#include "../Utils_Fill.h"

Functions
int	main ()

Detailed Description

Solves the harmonic equation.

$f''(x) + \lambda f(x) = 0$

as an eigenvalue problem for $\lambda$ over the unit domain with homogeneous boundary conditions for , returning any eigenvalue(s) with absolute value less than 10. The resulting eigenvalues should approach $m^2 \pi^2$ as $n \to \infty$ with $O(\Delta^2)$ corrections where $\Delta = 1/(n-1)$ and is the number of nodal points in the FD representation; a 2nd order central difference representation of

$f''(x_i) + \lambda f(x_i) = \frac{ f_{i-1} - 2f_i + f_{i+1} }{ \Delta^2 } + \lambda f_i + O(\Delta^2)$

is used. The matrix problem is solved for all eigenvalues within a specified distance of the origin of the comlplex plane by calling the LAPACK generalised eigenvalue routine.

Definition in file EVPHarmonic_lapack.cpp.

Function Documentation

◆ main()

int main ( )

Definition at line 25 of file EVPHarmonic_lapack.cpp.

{
  cout << "\n";
  cout << "=== EVP: Harmonic equation solved using LAPACK  =====\n";
  cout << "===  with a manually assembled matrix problem.\n";
  cout << "\n";
 
  cout.precision( 12 );
  cout << " Number of nodal points : Leading eigenvalue error : Total CPU time taken (ms) \n";
  bool failed = false;
  size_t N = 4;
  // a vector for the eigenvalues
  DenseVector<D_complex> lambdas;
  Timer timer;
  for ( int i = 2; i < 11; ++i )
  {
    N = ( size_t ) ( std::pow( 2., i ) );
    const double delta = 1. / ( N - 1 );
    const double delta2 = delta * delta;
    // matrix problem
    DenseMatrix<double> a( N, N, 0.0 );
    // Finite difference representation of f''(x)
    // here it's a a tri-diagonal system
    Utils_Fill::fill_band(a,-1, 1.0/delta2);
    Utils_Fill::fill_band(a, 0,-2.0/delta2);
    Utils_Fill::fill_band(a, 1, 1.0/delta2);
    // overwrite with boundary conditions at f(0) = f(1) = 0
    a( 0, 0 ) = 1.0;
    a( 0, 1 ) = 0.0;
    a( N - 1, N - 1 ) = 1.0;
    a( N - 1, N - 2 ) = 0.0;
    // not a generalised problem - but we'll apply that routine anyway
    // b is the RHS matrix, so it's -I
    DenseMatrix<double> b( N, N, 0.0 );
    Utils_Fill::fill_identity( b );
    b.scale( -1.0 );
    b( 0, 0 ) = 0.0;
    b( N - 1, N - 1 ) = 0.0;
    // a vector for the eigenvectors - although we won't use them
    DenseMatrix<D_complex> eigenvectors;
 
    DenseLinearEigenSystem<double> system( &a, &b );
    system.set_calc_eigenvectors( true );
 
    timer.start();
    try
    {
      system.eigensolve();
    }
    catch (const std::runtime_error &error )
    {
      cout << " \033[1;31;48m  * FAILED THROUGH EXCEPTION BEING RAISED \033[0m\n";
      return 1;
    }
    system.tag_eigenvalues_disc( + 1, 10. );
    lambdas = system.get_tagged_eigenvalues();
    eigenvectors = system.get_tagged_eigenvectors();
    cout << "    " << N << " : " << lambdas[ 0 ].real() - M_PI * M_PI
         << " : " << timer.get_time() << "\n";
    timer.stop();
  }
 
  const double tol = 1.e-4;
  if ( abs( lambdas[ 0 ].real() - M_PI * M_PI ) > tol )
    failed = true;
 
  if ( failed )
  {
    cout << "\033[1;31;48m  * FAILED \033[0m\n";
    cout << "    Final error = " << abs( lambdas[ 0 ].real() - M_PI * M_PI ) << "\n";
    return 1;
  }
 
  cout << "\033[1;32;48m  * PASSED \033[0m\n";
  return 0;
}

References CppNoddy::DenseLinearEigenSystem< _Type >::eigensolve(), Utils_Fill::fill_band(), Utils_Fill::fill_identity(), CppNoddy::DenseLinearEigenSystem< _Type >::get_tagged_eigenvalues(), CppNoddy::DenseLinearEigenSystem< _Type >::get_tagged_eigenvectors(), CppNoddy::Timer::get_time(), CppNoddy::DenseMatrix< _Type >::scale(), CppNoddy::LinearEigenSystem_base::set_calc_eigenvectors(), CppNoddy::Timer::start(), CppNoddy::Timer::stop(), and CppNoddy::DenseLinearEigenSystem< _Type >::tag_eigenvalues_disc().

Functions

Detailed Description

Function Documentation

◆ main()