Let's go through some application code a little at a time: The startup stuff is as follows:
#include "simulation.h"
#include "doublevec2.h"
int main(void)
{
FLOAT_POLY P,P2;
FLOAT_POLY Q;
float A[10] = {1,1, 2,2, 3,3, 4,4, 5,4};
float B[8] = {-1.0,1.0, 2.0,2.0, 3.0,3.0, 4.0,4.0};
Now let's build a polynomial:
Q.build_list(A,10); cout << "float Polynomial Q = " << endl << Q << endl;
This generates the output (we should fix our print method as the spacing in the output is less than perfect!):
float Polynomial Q = (1) x^9+(1) x^8+(2) x^7+(2) x^6+(3) x^5+(3) x^4+(4) x^3+(4) x^2+(5) x+(4) .
Now test the build_list method:
P.build_list(B,8);
cout << "P = " << endl << P << endl << endl;
This gives the output:
P = (-1) x^7+(1) x^6+(2) x^5+(2) x^4+(3) x^3+(3) x^2+(4) x+(4) .
Now the polynomial P2 was constructed using the empty or default constructor. So it starts with the head pointer NULL. Let's add some stuff to it:
P2.insert(1.2,6); P2.insert(2.4,4); P2.insert(3.5,1); P2.insert(4.5,0); cout << "Polynomial P2 = " << P2 << endl << endl;
generating output
Polynomial P2 = (1.2) x^6+(2.4) x^4+(3.5) x+(4.5) .
next, let's test polynomial addition, overloaded equal and so forth:
cout << "Test adding polynomials" << endl; cout << "P = " << P << endl; cout << "P2 = " << P2 << endl; cout << "Now add P + P2 and store in S" << endl; FLOAT_POLY S = P + P2; cout << "P + P2 = " << endl << S << endl << endl; cout << "Test adding polynomials" << endl; cout << "P = " << P << endl; cout << "Q = " << Q << endl; cout << "Now add P + Q and store in R" << endl; FLOAT_POLY R = P + Q; cout << "P + Q = " << endl << R << endl << endl;
We find
Test adding polynomials P = (-1) x^7+(1) x^6+(2) x^5+(2) x^4+(3) x^3+(3) x^2+(4) x+(4) . P2 = (1.2) x^6+(2.4) x^4+(3.5) x+(4.5) . Now add P + P2 and store in S P + P2 = (-1) x^7+(2.2) x^6+(2) x^5+(4.4) x^4+(3) x^3+(3) x^2+(7.5) x+(8.5) . Test adding polynomials P = (-1) x^7+(1) x^6+(2) x^5+(2) x^4+(3) x^3+(3) x^2+(4) x+(4) . Q = (1) x^9+(1) x^8+(2) x^7+(2) x^6+(3) x^5+(3) x^4+(4) x^3+(4) x^2+(5) x+(4) . Now add P + Q and store in R P + Q = (1) x^9+(1) x^8+(1) x^7+(3) x^6+(5) x^5+(5) x^4+(7) x^3+(7) x^2+(9) x+(8) .
Now check scalar multiplication:
R = (float)2.0*P2; cout << "P2 = " << P2 << endl; cout << "R = 2.0 * P2 = " << endl << R << endl << endl; R = P2*(float)2.0; cout << "P2 = " << P2 << endl; cout << "R = P2*2.0 = " << endl << R << endl << endl;
leading to the output:
P2 = (1.2) x^6+(2.4) x^4+(3.5) x+(4.5) .
R = 2.0 * P2 =
(2.4) x^6+(4.8) x^4+(7) x+(9) .
P2 = (1.2) x^6+(2.4) x^4+(3.5) x+(4.5) .
R = P2*2.0 =
(2.4) x^6+(4.8) x^4+(7) x+(9) .
Check the integration methods:
cout << "integral R =" << endl; cout << R.integral() << endl; cout << "definite integral from 1 to 2 of R = " << endl; cout << R.definite_integral(1.0,2.0) << endl;
with output
integral R = (0.342857) x^7+(0.96) x^5+(3.5) x^2+(9) x+(1) . definite integral from 1 to 2 of R = 92.8029
Now we will get more sophisticated and starting with a set of ten polynomials, {1, t, t2,..., t9}, we use Gram-Schmidt orthogonalization to compute ten new polynomials which are mutually orthogonal on the interval [- 1, 1] under the L2 norm. The code is
FLOAT_POLY *SET[10], BASIS[10],temp,temp2;
float norm;
int i,j;
for(i=0;i<10;++i){
SET[i] = new FLOAT_POLY(1.0,(unsigned int) i);
cout << "SET[" << i << "] = " << *SET[i] << endl;
cout << "SET[" << i << "](2.0) = " << (*SET[i])(2.0) << endl;
temp = SET[i]->integral();
cout << "SET[" << i << "].integral = " << SET[i]->integral() << endl;
cout << "SET[" << i << "].definite_integral(-1.0,1.0) = " <<
SET[i]->definite_integral(-1.0,1.0) << endl;
cout << endl;
}
This gives us the following polynomials:
SET[0] = (1) . SET[0](2.0) = 1 SET[0].integral = (1) x+(1) . SET[0].definite_integral(-1.0,1.0) = 2 SET[1] = (1) x. SET[1](2.0) = 2 SET[1].integral = (0.5) x^2+(1) . SET[1].definite_integral(-1.0,1.0) = 0 SET[2] = (1) x^2. SET[2](2.0) = 4 SET[2].integral = (0.333333) x^3+(1) . SET[2].definite_integral(-1.0,1.0) = 0.666667 SET[3] = (1) x^3. SET[3](2.0) = 8 SET[3].integral = (0.25) x^4+(1) . SET[3].definite_integral(-1.0,1.0) = 0 SET[4] = (1) x^4. SET[4](2.0) = 16 SET[4].integral = (0.2) x^5+(1) . SET[4].definite_integral(-1.0,1.0) = 0.4 SET[5] = (1) x^5. SET[5](2.0) = 32 SET[5].integral = (0.166667) x^6+(1) . SET[5].definite_integral(-1.0,1.0) = 0 SET[6] = (1) x^6. SET[6](2.0) = 64 SET[6].integral = (0.142857) x^7+(1) . SET[6].definite_integral(-1.0,1.0) = 0.285714 SET[7] = (1) x^7. SET[7](2.0) = 128 SET[7].integral = (0.125) x^8+(1) . SET[7].definite_integral(-1.0,1.0) = 0 SET[8] = (1) x^8. SET[8](2.0) = 256 SET[8].integral = (0.111111) x^9+(1) . SET[8].definite_integral(-1.0,1.0) = 0.222222 SET[9] = (1) x^9. SET[9](2.0) = 512 SET[9].integral = (0.1) x^10+(1) . SET[9].definite_integral(-1.0,1.0) = 0
The following code computes the new orthonormal basis from these original ten polynomials using the Gram-Schmidt orthogonalization procedure:
norm = SET[0]->inner_product(-1.0,1.0,*SET[0]);
norm = 1.0/sqrt(norm);
BASIS[0] = norm*(*SET[0]);
cout << "BASIS[0] = " << BASIS[0] << endl << endl;
for(i=1;i<10;++i){
temp = *SET[i];
for(j=0;j<i;++j){
float scale = SET[i]->inner_product(-1.0,1.0,BASIS[j]);
temp = temp - scale*BASIS[j];
}
norm = temp.inner_product(-1.0,1.0,temp);
norm = 1.0/sqrt(norm);
BASIS[i] = norm*temp;
BASIS[i].mergesort(); //sort the new polynomials
cout << "BASIS[" << i << "] = " << BASIS[i] << endl << endl;
}
generating the basis:
BASIS[0] = (0.707107) .
BASIS[1] = (1.22474) x+(-0) .
BASIS[2] = (2.37171) x^2+(-0) x+(-0.790569) .
BASIS[3] = (4.67707) x^3+(-0) x^2+(-2.80624) x+(0) .
BASIS[4] = (9.28076) x^4+(-0) x^3+(-7.95493) x^2+(0) x+(0.795494) .
BASIS[5] = (18.4684) x^5+(-0) x^4+(-20.5205) x^3+(0) x^2+(4.39724) x+(0) .
BASIS[6] = (36.8114) x^6+(-0) x^5+(-50.1971) x^4+(0) x^3+(16.7322) x^2
+(0) x+(-0.796765) .
BASIS[7] = (73.4244) x^7+(-0) x^6+(-118.608) x^5+(0) x^4+(53.9128) x^3
+(0) x^2+(-5.99027) x+(0) .
BASIS[8] = (146.286) x^8+(-0) x^7+(-273.126) x^6+(0) x^5+(157.619) x^4
+(0) x^3+(-28.6696) x^2+(0) x+(0.796856) .
BASIS[9] = (292.571) x^9+(-0) x^8+(-619.539) x^7+(0) x^6+(433.654) x^5
+(0) x^4+(-111.186) x^3+(0) x^2+(7.58003) x+(0) .
If we had inserted some additional prints into the for loop calculation, we would see all of the inner products used to calculate the basis functions:
BASIS[0] = (0.707107) . <(1) x., (0.707107) . > = 0 temp = (1) x. BASIS[1] = (1.22474) x. <(1) x^2., (0.707107) . > = 0.471405 temp = (1) x^2+ (-0.333333) . <(1) x^2., (1.22474) x. > = 0 temp = (1) x^2+ (-0.333333) . BASIS[2] = (2.37171) x^2+ (-0.790569) . <(1) x^3., (0.707107) . > = 0 temp = (1) x^3. <(1) x^3., (1.22474) x. > = 0.489898 temp = (1) x^3+ (-0.6) x. <(1) x^3., (2.37171) x^2+ (-0.790569) . > = 0 temp = (1) x^3+ (-0.6) x. BASIS[3] = (4.67707) x^3+ (-2.80624) x. ...
Now these basis functions should be mutually orthogonal and have dot product 1: hence if we compute the dot products
Basis[i].innner_product(-1.0,1.0,BASIS[j])
for all possible i and j indices, we should generate a 10 by 10 identity matrix.
The code to verify this is as follows:
DOUBLE_MATRIX T(10,10);
for(i=0;i<10;++i)
for(j=0;j<10;++j)
T[i][j] = BASIS[i].inner_product(-1.0,1.0,BASIS[j]);
cout << "Kronecker Delta Matrix = " << endl;
cout << T;
return 0;
}
with output (with a little eyestrain, you can see this is the 10 by 10 identity matrix as the entries on the order of 10-6 are roughly at machine zero for floats.
Kronecker Delta Matrix = 1 0 -5.96046e-08 0 7.15256e-07 0 -2.92063e-06 0 1.43051e-05 0 0 1 0 2.38419e-07 0 2.38419e-07 0 8.10623e-06 0 -2.24113e-05 -5.96046e-08 0 1 0 -5.36442e-07 0 -3.09944e-06 0 1.62721e-05 0 0 2.38419e-07 0 1 0 9.53674e-07 0 -4.76837e-06 0 -7.9155e-05 7.15256e-07 0 -5.36442e-07 0 0.999995 0 2.47359e-05 0 -0.0001266 0 0 2.38419e-07 0 9.53674e-07 0 1.00002 0 2.28882e-05 0 -2.09808e-05 -2.92063e-06 0 -3.09944e-06 0 3.99947e-05 0 1.00019 0 -0.00135899 0 0 8.10623e-06 0 -4.76837e-06 0 2.28882e-05 0 0.99996 0 -8.96454e-05 1.43051e-05 0 1.62721e-05 0 -0.0001266 0 -0.0017252 0 0.996796 0 0 -2.24113e-05 0 -7.9155e-05 0 -0.000143051 0 -0.000822067 0 0.992729