Numerical Methods for Physics - Program

Programs from the textbook "Numerical Methods for Physics" by Alejandro
Garcia.

Cpp/Matrix.h

@@ -0,0 +1,93 @@
+#include <assert.h>  // Defines the assert function.
+
+class Matrix {
+
+public:
+
+// Default Constructor. Creates a 1 by 1 matrix; sets value to zero. 
+Matrix () {
+  nRow_ = 1; nCol_ = 1;
+  data_ = new double [1];  // Allocate memory
+  set(0.0);                // Set value of data_[0] to 0.0
+}
+
+// Regular Constructor. Creates an nR by nC matrix; sets values to zero.
+// If number of columns is not specified, it is set to 1.
+Matrix(int nR, int nC = 1) {
+  assert(nR > 0 && nC > 0);    // Check that nC and nR both > 0.
+  nRow_ = nR; nCol_ = nC;
+  data_ = new double [nR*nC];  // Allocate memory
+  assert(data_ != 0);          // Check that memory was allocated
+  set(0.0);                    // Set values of data_[] to 0.0
+}
+
+// Copy Constructor.
+// Used when a copy of an object is produced 
+// (e.g., passing to a function by value)
+Matrix(const Matrix& mat) {
+  this->copy(mat);   // Call private copy function.
+}
+
+// Destructor. Called when a Matrix object goes out of scope.
+~Matrix() {
+  delete [] data_;   // Release allocated memory
+}
+
+// Assignment operator function.
+// Overloads the equal sign operator to work with
+// Matrix objects.
+Matrix& operator=(const Matrix& mat) {
+  if( this == &mat ) return *this;  // If two sides equal, do nothing.
+  delete [] data_;                  // Delete data on left hand side
+  this->copy(mat);                  // Copy right hand side to l.h.s.
+  return *this;
+}
+
+// Simple "get" functions. Return number of rows or columns.
+int nRow() const { return nRow_; }
+int nCol() const { return nCol_; }
+
+// Parenthesis operator function.
+// Allows access to values of Matrix via (i,j) pair.
+// Example: a(1,1) = 2*b(2,3); 
+// If column is unspecified, take as 1.
+double& operator() (int i, int j = 1) {
+  assert(i > 0 && i <= nRow_);          // Bounds checking for rows
+  assert(j > 0 && j <= nCol_);          // Bounds checking for columns
+  return data_[ nCol_*(i-1) + (j-1) ];  // Access appropriate value
+}
+
+// Parenthesis operator function (const version).
+const double& operator() (int i, int j = 1) const{
+  assert(i > 0 && i <= nRow_);          // Bounds checking for rows
+  assert(j > 0 && j <= nCol_);          // Bounds checking for columns
+  return data_[ nCol_*(i-1) + (j-1) ];  // Access appropriate value
+}
+
+// Set function. Sets all elements of a matrix to a given value.
+void set(double value) {
+  int i, iData = nRow_*nCol_;
+  for( i=0; i<iData; i++ )
+    data_[i] = value;
+}
+
+//*********************************************************************
+private:
+
+// Matrix data.
+int nRow_, nCol_;  // Number of rows, columns
+double* data_;     // Pointer used to allocate memory for data.
+
+// Private copy function.
+// Copies values from one Matrix object to another.
+void copy(const Matrix& mat) {
+  nRow_ = mat.nRow_;
+  nCol_ = mat.nCol_;
+  int i, iData = nRow_*nCol_;
+  data_ = new double [iData];
+  for(i = 0; i<iData; i++ )
+    data_[i] = mat.data_[i];
+}
+
+}; // Class Matrix
+

Cpp/NumMeth.h

@@ -0,0 +1,8 @@
+// General header file for C++ programs 
+// in "Numerical Methods for Physics" 
+
+#include <iostream.h>
+#include <fstream.h>
+#include <assert.h>  
+#include <math.h>
+#include "Matrix.h"

Cpp/SampList.h

@@ -0,0 +1,53 @@
+class SampList {
+
+public:
+
+// Class data (sorting lists)
+int ncell, nsamp;
+double *ave_n, *ave_ux, *ave_uy, *ave_uz, *ave_T;
+
+// Default Constructor. 
+SampList() {
+  initLists(1);
+}
+
+// Regular Constructor. 
+SampList(int ncell_in) {
+  initLists(ncell_in);
+}
+
+
+// Destructor. Called when a SampList object goes out of scope.
+~SampList() {
+  delete [] ave_n;   // Release allocated memory
+  delete [] ave_ux;
+  delete [] ave_uy;
+  delete [] ave_uz;
+  delete [] ave_T;
+}
+
+//*********************************************************
+
+private:
+
+// Initialization routine
+void initLists(int ncell_in) {
+  ncell = ncell_in;
+  nsamp = 0;
+  ave_n = new double [ncell+1];  // Allocate memory
+  ave_ux = new double [ncell+1];
+  ave_uy = new double [ncell+1];
+  ave_uz = new double [ncell+1];
+  ave_T = new double [ncell+1];
+  int i;
+  for( i=1; i<=ncell; i++ ) {
+	ave_n[i] = 0;
+	ave_ux[i] = 0;
+	ave_uy[i] = 0;
+	ave_uz[i] = 0;
+	ave_T[i] = 0;
+  }
+}
+
+}; // Class SampList
+

Cpp/SortList.h

@@ -0,0 +1,47 @@
+class SortList {
+
+public:
+
+// Class data (sorting lists)
+int ncell, npart, *cell_n, *index, *Xref;
+
+// Default Constructor. 
+SortList() {
+  initLists(1,1);
+}
+
+// Regular Constructor. 
+SortList(int ncell_in, int npart_in) {
+  initLists(ncell_in,npart_in);
+}
+
+
+// Destructor. Called when a SortList object goes out of scope.
+~SortList() {
+  delete [] cell_n;   // Release allocated memory
+  delete [] index;
+  delete [] Xref;
+}
+
+//*********************************************************
+
+private:
+
+// Initialization routine
+void initLists(int ncell_in, int npart_in) {
+  ncell = ncell_in;
+  npart = npart_in;
+  cell_n = new int [ncell+1];  // Allocate memory
+  index = new int [ncell+1];
+  Xref = new int [npart+1];
+  int i;
+  for( i=1; i<=ncell; i++ ) {
+	cell_n[i] = 0;
+	index[i] = 0;
+  }
+  for( i=1; i<=npart; i++ )
+	Xref[i] = 0;
+}
+
+}; // Class SortList
+

Cpp/advect.cpp

@@ -0,0 +1,109 @@
+// advect - Program to solve the advection equation 
+// using the various hyperbolic PDE schemes
+#include "NumMeth.h"
+
+void main() {
+
+  //* Select numerical parameters (time step, grid spacing, etc.).
+  cout << "Choose a numerical method: 1) FTCS, 2) Lax, 3) Lax-Wendroff : ";
+  int method; cin >> method;
+  cout << "Enter number of grid points: "; int N; cin >> N;
+  double L = 1.;     // System size
+  double h = L/N;    // Grid spacing
+  double c = 1;      // Wave speed
+  cout << "Time for wave to move one grid spacing is " << h/c << endl;
+  cout << "Enter time step: "; double tau; cin >> tau;
+  double coeff = -c*tau/(2.*h);    // Coefficient used by all schemes
+  double coefflw = 2*coeff*coeff;  // Coefficient used by L-W scheme
+  cout << "Wave circles system in " << L/(c*tau) << " steps" << endl;
+  cout << "Enter number of steps: "; int nStep; cin >> nStep;
+
+  //* Set initial and boundary conditions.
+  const double pi = 3.141592654;
+  double sigma = 0.1;              // Width of the Gaussian pulse
+  double k_wave = pi/sigma;        // Wave number of the cosine
+  Matrix x(N), a(N), a_new(N);
+  int i,j;
+  for( i=1; i<=N; i++ ) {
+    x(i) = (i-0.5)*h - L/2;  // Coordinates of grid points
+    // Initial condition is a Gaussian-cosine pulse
+    a(i) = cos(k_wave*x(i)) * exp(-x(i)*x(i)/(2*sigma*sigma)); 
+  }
+  // Use periodic boundary conditions
+  int *ip, *im;  ip = new int [N+1];  im = new int [N+1];
+  for( i=2; i<N; i++ ) {
+    ip[i] = i+1;    // ip[i] = i+1 with periodic b.c.
+    im[i] = i-1;    // im[i] = i-1 with periodic b.c.
+  }
+  ip[1] = 2;  ip[N] = 1;
+  im[1] = N;  im[N] = N-1;
+
+  //* Initialize plotting variables.
+  int iplot = 1;          // Plot counter
+  int nplots = 50;        // Desired number of plots
+  double plotStep = ((double)nStep)/nplots; 
+  Matrix aplot(N,nplots+1), tplot(nplots+1);
+  tplot(1) = 0;       // Record the initial time (t=0)
+  for( i=1; i<=N; i++ )
+    aplot(i,1) = a(i);  // Record the initial state
+
+  //* Loop over desired number of steps.
+  int iStep;
+  for( iStep=1; iStep<=nStep; iStep++ ) {  
+
+    //* Compute new values of wave amplitude using FTCS, 
+    //  Lax or Lax-Wendroff method.
+    if( method == 1 )      ////// FTCS method //////
+      for( i=1; i<=N; i++ )
+        a_new(i) = a(i) + coeff*( a(ip[i])-a(im[i]) );  
+    else if( method == 2 )  ////// Lax method //////
+      for( i=1; i<=N; i++ )
+        a_new(i) = 0.5*( a(ip[i])+a(im[i]) ) + 
+                   coeff*( a(ip[i])-a(im[i]) );   
+    else                   ////// Lax-Wendroff method //////
+      for( i=1; i<=N; i++ )
+        a_new(i) = a(i) + coeff*( a(ip[i])-a(im[i]) ) +
+                   coefflw*( a(ip[i])+a(im[i])-2*a(i) ); 
+
+    a = a_new;	 // Reset with new amplitude values
+
+    //* Periodically record a(t) for plotting.
+    if( fmod((double)iStep,plotStep) < 1 ) {  
+      iplot++;
+      tplot(iplot) = tau*iStep;
+      for( i=1; i<=N; i++ )    
+        aplot(i,iplot) = a(i);       // Record a(i) for ploting
+      cout << iStep << " out of " << nStep << " steps completed" << endl;
+    }
+  }
+  nplots = iplot;   // Actual number of plots recorded
+
+  //* Print out the plotting variables: x, a, tplot, aplot
+  ofstream xOut("x.txt"), aOut("a.txt"), 
+	      tplotOut("tplot.txt"), aplotOut("aplot.txt");
+  for( i=1; i<=N; i++ ) {
+    xOut << x(i) << endl;
+    aOut << a(i) << endl;
+    for( j=1; j<nplots; j++ )
+      aplotOut << aplot(i,j) << ", ";
+    aplotOut << aplot(i,nplots) << endl;
+  }
+  for( i=1; i<=nplots; i++ )
+    tplotOut << tplot(i) << endl;
+
+  delete [] ip, im;  // Release allocated memory
+}
+/***** To plot in MATLAB; use the script below ********************
+%* Plot the initial and final states.
+load x.txt; load a.txt; load tplot.txt; load aplot.txt;
+figure(1); clf;  % Clear figure 1 window and bring forward
+plot(x,aplot(:,1),'-',x,a,'--');
+legend('Initial','Final');
+xlabel('x');  ylabel('a(x,t)');
+pause(1);    % Pause 1 second between plots
+%* Plot the wave amplitude versus position and time
+figure(2); clf;  % Clear figure 2 window and bring forward
+mesh(tplot,x,aplot);
+ylabel('Position');  xlabel('Time'); zlabel('Amplitude');
+view([-70 50]);  % Better view from this angle
+******************************************************************/