Numerical Integration of Acoustic Energy

Corey Mutnik
Physics 305 Computational Physics
2/3/2015

Introduction
Integral Estimations
Code
Graphs
Analysis
Conclusion
References
Acknowledgments
Home

1. Introduction

    Many integrals are to difficult to solve out analytically.  One such example is the chirp function:

                                                                                                                                                     (1)

Equation 1 has no known real anti-derivative.  The best way to solve the integral of such a function is to break the area under the generated curve into small portions that can then be summed up.  The closer a desired answer is to the actual result depends on how the segments, that make up the area beneath the curve, are chosen.  Ideally one would use segments of height f(t) and width "dt" and sum the area of the corresponding rectangles as "dt" goes to zero.  In order to do this am infinite number of segments would need to be summed.  Since this can not be done we use approximations.  Calculating and summing the area of each successive segment is a tedious and painstaking task.  Luckily computers are able to sum much larger numbers of iterations, at speeds factorially faster than had they been done by hand. 

    Computational analysis is extraordinarily helpful in various fields.  This lab will focus on its usefulness related to numerical integration of acoustic energy.
   

    The chirp function with an exponential damping term can be seen in the graph above.  It is apparent that the function is symmetric about t = 0, and rapidly oscillates with an amplitude of 2.  The damping term acts as a limiting agent, causing the function to quickly converge to zero just before t = +/- 6 (as shown in the plot above).  This allows us to measure the acoustic energy of the chirp function by integrating from -6 to 6, rather than using one that extends from negative infinity to infinity.

    Acoustic wave energy is a very complex thing to analyze.  "For simplicity here we will consider plane acoustic waves with no absorption."[1]  The propagation of sound waves occurs as longitudinal pressure waves that have an intensity:
                                                                                                                                       (2)

By integrating the third and fourth terms in equation 2, with respect to time, we obtain the energy transferred per unit area.

    When speaking of sound one often hears the term decibels (dB).  Decibels are not a unit of measurement on their own, they are relative.  That is to say that dB need a reference point to be compared to.  For acoustic intensity and pressure we have the relation:

                                                                                                                         (3)
Where Iref = 1 picoWatt/m^2 and  Pref = 20 microPascals.


2. Integral Estimations

    This lab implements two methods of estimating integrals by summing up segmented portions from the area between the curve and the axis. The Trapezoidal method and Simpson's method were chosen for their simplicity.  Each method was used to estimate the integral value of the chirp function with a damping term attached:

                                                                                                                                                (4)

The inability of the chirp function to be solved analytically makes it a prime candidate for computational analysis.  Here we allow the computer to sum extremely high numbers of intervals, in an attempt to achieve a more accurate result.  These methods allow us to estimate the value of the chirp function with extremely high precision.

Trapezoidal Method [3]

                                                    \int_{a}^{b} f(x)\, dx \approx \frac{h}{2} \sum_{k=1}^{N} \left( f(x_{k+1}) + f(x_{k}) \right)                                                               (5)
   
    Where h is the width of each interval.  The Trapezoidal method breaks each segment up into a trapezoid then sums them together to estimate the integral.  The trapezoidal method converges quickly for periodic functions, such as cosine.

Simpson's Method [4]

                                     \int_a^b f(x) \, dx\approx\tfrac{h}{3}\bigg[f(x_0)+2\sum_{j=1}^{n/2-1}f(x_{2j})+ 4\sum_{j=1}^{n/2}f(x_{2j-1})+f(x_n) \bigg],                               (6)

    Where h is the width of each interval.  Simpson's method applies a local quadratic fit to the function at each segment.  It then sums all the segments together to estimate the integral.


3. Code


Program 1
using namespace std;      
#include <iostream>      
#include <iomanip>   // this is required to use "setprecision" below
#include <fstream>
#define _USE_MATH_DEFINES 
#include <cmath>          

#define max_in  3001                  // max number of intervals
#define vmin 0.0                        // ranges of integration
#define vmax 6.0           
//#define ME 2.7182818284590452354E0      // Euler's number
// Since we used the math defines above, we can just use M_E for Euler's number...
#define ME M_E

float trapez  (int no, float min, float max);    // trapezoid rule
float simpson (int no, float min, float max);    // Simpson's rule
float f(float);                    //  contains function to integrate

main()
{         
   float result; // declare this within scope of main
   ofstream myfile;  // declare a new instance of an output stream, call it "myfile"
  
 
   myfile.open("integ.dat");  // open is now a method within ofstream class instance outfile
  
   for (int i=3; i<=max_in; i+=2){       // Simpson's rule requires odd number of intervals
      result = trapez(i, vmin, vmax);
      myfile << i <<"\t"<< setprecision(9) << result <<"\t";  // "\t" here gives a tab
      result = simpson(i, vmin, vmax);
      myfile << setprecision(9) << result << endl;
   }
   cout << "data stored in integ.dat" << endl;
   myfile.close();  // close the file, disable any more writing to it, button it up.
}
/*------------------------end of main program-------------------------*/

// the function we want to integrate
float f (float x)
{
return cos(pow(x,3.))*exp(-pow(x/4.,8.));  // this function for chirp
//return (exp(-x));    // this is exponential
}

/* Integration using trapezoid rule */
float trapez (int no, float min, float max)
{           
   float interval, sum=0., x;         

   interval = ((max-min) / (no-1));  // this is equivalent to "dt" or "dx" in integral
  
   sum += 0.5 * f(min) * interval; // add the first point
   for (int n=2; n<no; n++)               // sum the midpoints
   {
      x = min + interval * (n-1);     
      sum += f(x)*interval;
   }
   sum += 0.5 *f(max) * interval;    // add the endpoint

   return (sum);
}     

/* Integration using Simpson's rule */
float simpson (int no, float min, float max)
{                   
   float interval, sum=0., x;
   interval = ((max -min) /(no-1));
  
   for (int n=2; n<no; n+=2)                // loop for odd points 
   {
       x = min + interval * (n-1);
       sum += 4 * f(x);
   }
   for (int n=3; n<no; n+=2)                // loop for even points 
   {
      x = min + interval * (n-1);
      sum += 2 * f(x);
   }  
   sum +=  f(min) + f(max);         // add first and last value         
   sum *= interval/3.;                // then multilpy by interval
  
   return (sum);
}

    Program 1 was used to output values, corresponding to the chirp function, into a data file.  From their the information was plotted in gnuplot.  This code was also used to generate the plot for the dampened chirp function and the dampened chirp function with cos (x^3) replaced by cos(x^5).  Commenting out and in the proper functions allowed this code to be used in generating values for the exponential too.  Program 1 calculated integrals from 0 to 6 using both the Trapezoidal and Simpson's methods; with the number of terms being 3001 (max_in).


Program 2
// Integrate acoustic energy

using namespace std;      
#include <iostream>      
#include <iomanip>   // this is required to use "setprecision" below
#include <fstream>
#define _USE_MATH_DEFINES 
#include <cmath>          

#define max_in  600000                  // sets max number of intervals
#define vmin -6.0                        // sets lower limit of integration
#define vmax 6.0            // sets upper limit of integration
#define ME M_E                // Euler's number

# define RHO 1.2            // RHO is the density of air in kg/m^(3)
# define c 340.0            // c is the speed of sound in m/s
# define AMP 20.0            // A is the amplitude of the pressure pulse in Pa

double trapez (int no, double min, double max);        // trapezoid rule
double psq(double);                    // pressure squared
ofstream myfile;                      // declare a new instance of an output stream, call it "myfile"

main()
{
    double result;
    myfile.open("myinteg4.dat");              // open is now a method within ofstream class instance outfile
    result = trapez(max_in, vmin, vmax);
    cout << "Answer= " << setprecision(12) << result << endl;
    cerr << "data stored in myinteg4.dat" << endl;
  myfile.close();                    // close the file, disable any more writing to it, button it up
}


/*------------------------end of main program-------------------------*/


// the function we want to integrate
double psq (double t)

{
return pow(cos(pow(2*M_PI*t,3.))*exp(-pow(t/4.,8.))*AMP,2.);  // this function for chirp where 2pi is omega in rad/sec and t is time
}


/* Integration using trapezoid rule */
double trapez (int no, double min, double max)
{           
   double interval, sum=0.0, t;         

   interval = ((max-min) / (no-1.0));  // interval is same as "dt" in integral, max-min is change in time (delta t), no-1. is number of samples
   sum += 0.5 * psq(min) * 1.0/(RHO*c) * interval;     // add the first point
   for (int n=2; n<no; n++)                       // sum the midpoints
   {
      t = min + interval * (n-1.0);     
      sum += psq(t) * 1.0/(RHO*c) * interval;
      myfile << t << "\t" << setprecision(12) << sum << endl;        // "\t" here gives a tab
   }
   sum += 0.5 * psq(max) * 1.0/(RHO*c) * interval;            // add the endpoint

   return (sum);
} 

    Program 2 used the Trapezoidal method to calculate the integral of the function from -6 to 6, using 600,000 terms.  It has already been modified to print out the running sum for the trapez function as the sum is accumulating: shown by figure 4.


4.  Graphs



Figure 1: Error when Calculating Exponential Function



Figure 2: Comparison of methods for integrating damped chirp function; using t^3                                       

Figure 3: Comparison of methods for integrating damped chirp function; using t^5  



Figure 4: Shows the cumulative value vs time

5. Analysis

    Figure 1 shows how the Trapezoidal and Simpson methods the integral of the exponential function: e^(-x).  Both methods quickly fall to zero, as the number of terms increases.  This is expected since figure 1 represents the error each method holds with respect to the number of terms being added.

    Figures 2 and 3 show how the Trapezoidal and Simpson's methods of integrating the dampened chirp function compare.  Figure 2 uses the chirp function including the exponential damping term.  Figure 3 uses the same integrand with t^3 replaced by t^5.  The largest difference in these figures is how many intervals they take before converging.  Both figure 2 and 3 appear noisy at first since each successive term added in the integration causes a drastic variation in the overall sum (subtractive cancellation - alternating signs).  As the number of intervals, n, increases each integral converges to some value.  In both figures the Trapezoidal method converges to a stable answer quicker than the Simpson method.  This is in direct conflict with the prediction made in "Computational Physics" by Landau and Paez [1]; which claims the Simpson's method should be more accurate.  Since the Trapezoidal method converges faster it is the one chosen to calculate the total energy per unit area.

     The total energy per unit area delivered by the entire pulse is 3.48236979162 Joules per meter squared.  Figure 4 shows how the energy per unit area changes over time.  The acoustic intensity (time rate of change of energy per unit area) is the derivative of this function.  Figure 4 shows that the energy per unit area has an inflection point about t=0, where it changes from an increasing to decreasing function.  Therefore we can conclude the intensity has a maximum at t=0.  The time rate of change of delivered energy is at a maximum halfway through the interval in which the pulse occurs.

    Using equation 3 we are able to calculate the acoustic dB during the pulse.  First we divide our total energy per unit area by the time interval and substitute it into our equation as our Intensity, I.  Remembering to convert Iref  = 1 picoWatt = 1e-12 J/s, allows our units to cancel.  Giving us a result of dB = 114.67 during the pulse.


6. Conclusions

  Computational methods of numerical analysis allow for much high accuracy in a fraction of the time.  The program written to estimate the value of a non-analytic function gave extremely precise results almost instantaneously.  This could not have been done by hand, or using any other method.  Thus proving computational numerical analysis to be an essential tool in the fields of physics and mathematics.

 

7. References


[1] R. H. Landau & M. J. Paez, Computational Physics, Problem Solving with Computers,  (Wiley: New York)                1997.

[2] Equations generated using: http://www.codecogs.com/latex/eqneditor.php

[3] http://www.mathwords.com/t/trapezoid_rule.htm

[4] http://www.mathwords.com/s/simpsons_rule.htm
 
[5]  "Physical Audio Signal Processing'', by Julius O. Smith III, W3K Publishing, 2010, ISBN                
                978-0-9745607-2-4.Copyright © 2014-10-15 by Julius O. Smith III


8. Acknowledgments 


1. Alexander, Ryan

2. Gorham, Peter

Home - Lab 1 - Lab 2 - Lab 3