Answer to problem 6.1

Below is the program and the output of the Gauss-Laguerre approximation for the integral of Plank's black body radiation equation.  This program verified that as n grows, the integral from zero to infinity of x^3/(e^x-1) = (p^4)/15.  Given the total power radiated from a black body for all frequencies.

we can substitute v=kTx/h to achieve an integral equivalent to the one evaluated times a factor 2ph*(kT/h)^3*(kT/h).  Plugging our approximation into the equation we achieve

Rt = (2ph/c^2)*(kT/h)^3*(kT/h)*(p^4/15) = (2(k^4)(p^5) / 15c^2) T^4 = sT^4

gausslaguerre.cpp

#include <functional>
#include <iostream>
#include <math.h>
#include <vector>
using namespace std;

const double pi = acos(-1);

double xm2[2] = {5.857864376269050e-1, 3.414213562373095};
double wm2[2] = {8.535533905932738e-1,1.464466094067262e-1};
double xm4[4] = {3.225476896193923e-1,1.745761101158347,  4.536620296921128,  9.395070912301133};
double wm4[4] = {6.031541043416336e-1,  3.574186924377997e-1,  3.888790851500538e-2,  5.392947055613275e-4 };
double xm6[6] = {2.228466041792607e-1 , 1.188932101672623    , 2.992736326059314     , 5.775143569104511   , 9.837467418382590   , 1.598287398060170e-1};
double wm6[6] = {4.589646739499636e-1  ,  4.170008307721210e-1   , 1.133733820740450e-1  , 1.039919745314907e-2   , 2.610172028149321e-4   , 8.985479064296212e-7};
double xmwm8[8][3] = 
{ 
	0.170279632305, 0.369188589342, 0.437723410493 ,	
	0.903701776799, 0.418786780814, 1.03386934767 ,
	2.25108662987, 0.175794986637 ,1.66970976566 ,
	4.26670017029, 0.0333434922612 ,2.37692470176, 
	7.04590540239, 0.00279453623523 ,3.20854091335 ,
	10.7585160102, 9.07650877338E-005, 4.26857551083 ,
	15.7406786413, 8.48574671626E-007 ,5.81808336867 ,
	22.8631317369, 1.04800117487E-009 ,8.90622621529 
};

class Function : public unary_function<double,double>
{
public:
	Function() {};
	virtual result_type operator()(argument_type x){return x;};
};

class BlackBody : public Function
{
public:
	BlackBody() {};
	double operator()(double x)
	{
		return x*x*x/(1-exp(-x));
	}
};

class Approximation
{
public:
	Approximation(double lowerBound, double upperBound, Function* function)
		:a(lowerBound), b(upperBound), f(function)
	{};

	virtual double calculate(double precision=0, int maxTries = 3000)
	{
		cerr << "Default approximation not implemented" << endl;
		return -1;
	}
protected:
	double a;
	double b;
	Function* f;
	//TBD: set precision & maxTries?
};

class GaussLaguerre : public Approximation
{
public:
	GaussLaguerre(int numberOfPoints, Function* f)	:Approximation(1, numberOfPoints, f), weights(8), abscissas(8)
	{
		int i;
		switch(numberOfPoints)
		{
		case 2:
			for (i=0;i<2;i++)
			{
				weights[i] = wm2[i];
				abscissas[i] = xm2[i];
			}
			break;
		case 4:
			for (i=0;i<4;i++)
			{
				weights[i] = wm4[i];
				abscissas[i] = xm4[i];
			}
			break;
		case 6:
			for (i=0;i<6;i++)
			{
				weights[i] = wm6[i];
				abscissas[i] = xm6[i];
			}
			break;
		case 8:
			for (i=0;i<8;i++)
			{
				weights[i] = xmwm8[i][1];
				abscissas[i] = xmwm8[i][0];
				//cout << abscissas[i] << "\t" << weights[i] << endl;
			}
			break;
		};	
	}

	virtual double calculate()
	{
		double sum=0;
		for (int i =a;i<=b;i++)
			sum += (*f)(abscissas[i-1])*weights[i-1];
		return sum;
	}

protected:
	vector<double> weights;
	vector<double> abscissas;

};

int main(int argc, char* argv)
{
	cout.precision(15);
	cout << "Here are the Gauss Laguerre Approximations of Plank's Black Body equation where pi = " << pi << endl;
	cout << "N\tvalue\terror" << endl;
	BlackBody bb;
	for (int i =2;i<=8;i+=2)
	{
		GaussLaguerre gl(i,&bb);
		cout << i << "\t" << gl.calculate() << "\t" << fabs(gl.calculate() - pi*pi*pi*pi/15) <<endl;; 
	}
return 0;
}

laguerre.txt

Here are the Gauss Laguerre Approximations of Plank's Black Body equation where pi = 3.14159265358979
N value                         error
2 6.41372746951758 0.080211932749247
4 6.49453563980263 0.000596237535802402
6 6.49027279227319 0.00366660999364132
8 6.49393566536404 3.73690278987482e-006