any suggestions to improve code?

Question

sexyzebra19 0 Junior Poster in Training

14 Years Ago

Could someone please have a look at my code and tell me anything that can be improved? Anything at all - layout, using headers etc?

Thanks in advance!

#include<iostream>
#include<vector>
#include<cmath>
#include<iomanip>

using namespace std;

typedef double(*fun)(double);

//Declare functions
double gauss(fun f, double &a, double &b, int &n, double &t, double &g1, double &g2, double &g3, double &g4, double &g5, double &h1, double &h2, double &h3, double &h4, double &h5, double &answer, vector<double> &abscisass, vector<double> &weights);
double equation(double t);
double err(double &answer, double &exactanswer, double &error);


//=============================================================================================================


int main()
{
    double a = 0.0;
    double b = 6.0;
    double g1 = 0.0;
    double g2 = 0.0;
    double g3 = 0.0;
    double g4 = 0.0;
    double g5 = 0.0;
    double h1 = 0.0;
    double h2 = 0.0;
    double h3 = 0.0;
    double h4 = 0.0;
    double answer = 0.0;
    double exactanswer = 76.750557;
    double error = 0.0;
    int n = 0;
    double t = 0;
    vector<double> abscisass(30,0);
    vector<double> weights(30,0);


    // n = 1
    abscisass[0] = 0.0;
    weights  [0] = 2.000000000000000;
    // n = 2
    abscisass[1] = 0.577350269189626;
    weights  [1] = 1.000000000000000;
    // n = 3
    abscisass[2] = 0.0;
    weights  [2] = 0.888888888888889;
    abscisass[3] = 0.774596669241483;
    weights  [3] = 0.555555555555556;
    // n = 4
    abscisass[4] = 0.339981043584856;
    weights  [4] = 0.652145154862546;
    abscisass[5] = 0.861136311594053;
    weights  [5] = 0.347854845137454;
    // n = 5
    abscisass[6] = 0.0;
    weights  [6] = 0.568888888888889;
    abscisass[7] = 0.538469310105683;
    weights  [7] = 0.478628670499366;
    abscisass[8] = 0.906179845938664;
    weights  [8] = 0.236926885056189;

  gauss(equation,a,b,n,t,g1,g2,g3,g4,g5,g6,g7,g8,g9,g10,h1,h2,h3,h4,h5,h6,h7,h8,h9,h10,answer,abscisass,weights);

    err(answer,exactanswer,error);

    cout << "the answer is " << setprecision(20) << answer << endl;

    cout << "the error is " << setprecision(6) << error << endl;

    return 0;
}

double gauss(fun f, double &a, double &b, int &n, double &t, double &g1, double &g2, double &g3, double &g4, double &g5, double &h1, double &h2, double &h3, double &h4, double &h5, double &answer, vector<double> &abscisass, vector<double> &weights)
{
    cout << "\n********** GAUSS QUADRATURE **********"<< endl;
    cout << "For the equation f(t) = te^2t" << endl;
    cout << "Enter the quadrature order, n between 1 and 10 you wish to use: " << endl;
    cin >> n;

    if(n==1)
    {
        t = abscisass[0];
        g1 = (((b - a)/2)*t) + (b + a)/2;      //coordinate transformation
        h1 = (weights[0])*((b - a)/2)*(f(g1));              // substitution into original equation
        answer = h1;
    }

    if(n==2)
    {
        t = abscisass[1];
        g1 = (((b - a)/2)*t) + (b + a)/2;      //coordinate transformation
        h1 = (weights[1])*((b - a)/2)*(f(g1));              // substitution into original equation
        t = -1*abscisass[1];
        g2 = (((b - a)/2)*t) + (b + a)/2;
        h2 = (weights[1])*((b - a)/2)*(f(g2));
        answer = h1 + h2;
    }

    else if(n==3)
    {
        t = abscisass[2];
        g1 = (((b - a)/2)*t) + (b + a)/2;      //coordinate transformation
        h1 = (weights[2])*((b - a)/2)*(f(g1));              // substitution into original equation
        t = abscisass[3];
        g2 = (((b - a)/2)*t) + (b + a)/2;
        h2 = (weights[3])*((b - a)/2)*(f(g2));
        t = -1*abscisass[3];
        g3 = (((b - a)/2)*t) + (b + a)/2;
        h3 = (weights[3])*((b - a)/2)*(f(g3));
        answer = h1 + h2 + h3;
    }
    else if(n==4)
    {
        t = abscisass[4];
        g1 = (((b - a)/2)*t) + (b + a)/2;      //coordinate transformation
        h1 = (weights[4])*((b - a)/2)*(f(g1));              // substitution into original equation
        t = -1*abscisass[4];
        g2 = (((b - a)/2)*t) + (b + a)/2;
        h2 = (weights[4])*((b - a)/2)*(f(g2));
        t = abscisass[5];
        g3 = (((b - a)/2)*t) + (b + a)/2;
        h3 = (weights[5])*((b - a)/2)*(f(g3));
        t = -1*abscisass[5];
        g4 = (((b - a)/2)*t) + (b + a)/2;
        h4 = (weights[5])*((b - a)/2)*(f(g4));
        answer = h1 + h2 + h3 + h4;
    }
    else if(n==5)
    {
        t = abscisass[6];
        g1 = (((b - a)/2)*t) + (b + a)/2;      //coordinate transformation
        h1 = (weights[6])*((b - a)/2)*(f(g1));    // substitution into original equation
        t = abscisass[7];
        g2 = (((b - a)/2)*t) + (b + a)/2;
        h2 = (weights[7])*((b - a)/2)*(f(g2));
        t = -1*abscisass[7];
        g3 = (((b - a)/2)*t) + (b + a)/2;
        h3 = (weights[7])*((b - a)/2)*(f(g3));
        t = abscisass[8];
        g4 = (((b - a)/2)*t) + (b + a)/2;
        h4 = (weights[8])*((b - a)/2)*(f(g4));
        t = -1*abscisass[8];
        g5 = (((b - a)/2)*t) + (b + a)/2;
        h5 = (weights[8])*((b - a)/2)*(f(g5));
        answer = h1 + h2 + h3 + h4 + h5;
    }

double equation (double t)
{
    return log(t) +4*t;
}

double err(double &answer, double &exactanswer, double &error)
{
    error = fabs(((answer - exactanswer)/exactanswer)*100);
}

c++

3 Contributors
11 Replies
131 Views
1 Day Discussion Span
Latest Post 14 Years Ago Latest Post by jonsca

jonsca 1,059 Quantitative Phrenologist

14 Years Ago

You shouldn't repeat all those calculations over and over in each branch of the if statements.
Do the common calculations, have an if statement that changes any values that are uncommon to the set, then proceed with the next set of common calculations for the group, etc.

Also, you should pass in arrays, structs, or vectors of your parameters to your gauss function to clean up the sheer volume of values. It's just a lot to have to put in there and once you get beyond 4-5 parameters it starts to make your head spin.

Other than that, seems to be ok but I didn't put it under a microscope. Keep plugging away at those numerical methods! :)

Dave Sinkula 2,398 long time no c

14 Years Ago

Could someone please have a look at my code and tell me anything that can be improved? Anything at all - layout, using headers etc?

Fix the errors?

main.cpp
main.cpp: In function `int main()':
main.cpp:65: error: `g6' was not declared in this scope
main.cpp:65: error: `g7' was not declared in this scope
main.cpp:65: error: `g8' was not declared in this scope
main.cpp:65: error: `g9' was not declared in this scope
main.cpp:65: error: `g10' was not declared in this scope
main.cpp:65: error: `h5' was not declared in this scope
main.cpp:65: error: `h6' was not declared in this scope
main.cpp:65: error: `h7' was not declared in this scope
main.cpp:65: error: `h8' was not declared in this scope
main.cpp:65: error: `h9' was not declared in this scope
main.cpp:65: error: `h10' was not declared in this scope
main.cpp:65: warning: unused variable 'g6'
main.cpp:65: warning: unused variable 'g7'
main.cpp:65: warning: unused variable 'g8'
main.cpp:65: warning: unused variable 'g9'
main.cpp:65: warning: unused variable 'g10'
main.cpp:65: warning: unused variable 'h5'
main.cpp:65: warning: unused variable 'h6'
main.cpp:65: warning: unused variable 'h7'
main.cpp:65: warning: unused variable 'h8'
main.cpp:65: warning: unused variable 'h9'
main.cpp:65: warning: unused variable 'h10'
main.cpp: In function `double gauss(double (*)(double), double&, double&, int&, double&, double&, double&, double&, double&, double&, double&, double&, double&, double&, double&, double&, std::vector<double, std::allocator<double> >&, std::vector<double, std::allocator<double> >&)':
main.cpp:152: error: a function-definition is not allowed here before '{' token
main.cpp:157: error: a function-definition is not allowed here before '{' token
main.cpp:159: error: expected `}' at end of input
*** Errors occurred during this build ***

jonsca commented: Point taken! :D +2

Dave Sinkula 2,398 long time no c

14 Years Ago

First, I guess, I'm not terribly fond of user-interactive programs -- the ones that pause and ask you to enter a value, etc. I think of your gauss() function as more a a math thing: given some input(s) it will produce some output(s). I find it much easier for testing to move any user-prompting out of this function (if to include it at all). But that's maybe just me.

I might have something that looks more like this:

#include <iostream>
#include <cmath>
#include <iomanip>
using namespace std;

double equation(double t)
{
   return log(t) + 4 * t;
}

double err(double answer, double exactanswer)
{
   return fabs(((answer - exactanswer) / exactanswer) * 100);
}

typedef double(*fun)(double);

double gauss(fun f, double a, double b, int n);

int main()
{
   double exactanswer = 76.750557;
   for ( int n = 1; n <= 10; ++n )
   {
      double answer = gauss(equation, 0, 6, n);
      cout << "n = " << setw(2) << n 
           << ": answer = " << setprecision(20) << answer 
           << ", error = "  << setprecision(6)  << err(answer, exactanswer) 
           << endl;
   }
   return 0;
}

Your original equation() function was a nice example: it takes a double and returns a double; no references, just simple. But then your err() function passed the result as a reference (even though it had a return value, which was unused/discarded later): I'd prefer to pass the two parameters and return the result. It seems more "natural" to me.

But it would seem that you have some predisposition to passing reference parameters, as noted in your version(s) of gauss(). Why? You're not changing these values in the called function for use in the caller so much as declaring temporaries for the function being called. This seems rather strange to me.

Why not just put the temporary variables in the function that uses them? Then you don't need to have such a huge function signature either. If you need more temporaries in gauss(), then put them there.

And the same thing for the constants. Why declare the constants that are used only by gauss() in some other function and then pass them to gauss(). Why not just put them in gauss(), since that is where they are used, and be done with it. I find it preferable to keep variables as local as possible.

If you get right down to it, then, your gauss() function doesn't really need very many parameters: the function to test, the a and b values, and the n value for the calculation. I'll get to that shortly.

In programming, you should generally develop the ability to spot patterns. Where patterns present themselves, we often take advantage of that in code that helps avoid repetition. For example,

else if(n==3)
    {
        t = abscisass[2];
        g = (((b - a)/2)*t) + (b + a)/2;      //coordinate transformation
        h = (weights[2])*((b - a)/2)*(f(g));              // substitution into original equation
        t = abscisass[3];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[3])*((b - a)/2)*(f(g));
        t = -1*abscisass[3];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[3])*((b - a)/2)*(f(g));
        answer = h;
    }

You've got (b-a)/2 and (b+a)/2 in a number of spots. Maybe it gets optimized by the compiler, but using a temporary makes it easier to follow as well as avoiding repeated recalculation.

double x = (b + a) / 2;
   double y = (b - 2) / 2;
   // ...
    else if(n==3)
    {
        t = abscisass[2];
        g =  y * t + x;      //coordinate transformation
        h = weights[2] * y * f(g);              // substitution into original equation
        t = abscisass[3];
        g = y * t + x;
        h = h + weights[3]  * y * f(g);
        t = -1*abscisass[3];
        g = y * t + x;
        h = h + weights[3] * y * f(g);
        answer = h;
    }

To me, this is a little easier to read and follow.

Another pattern: You calculate t, g, h. Then t, g, h. Then t, g, h. In C++ we have the inline function option (as opposed to a multiline macro in C) that would seem well suited to simplify some of this. Making a little helper function, much like using a temporary variable, can encapsulate this pattern in code.

double helper(fun f, double t, double w, double x, double y)
{
   double g = y * t + x;
   return w * y * f(g);
}

// ...

   case 3:
      h += helper(f,  table[2].abscisass, table[2].weights, x, y);
      h += helper(f,  table[3].abscisass, table[3].weights, x, y);
      h += helper(f, -table[3].abscisass, table[3].weights, x, y);
      break;

There seems to be something of a pattern to the n cases, but I didn't care to venture much further. The gauss() function, and its helper(), might then end up like this:

double helper(fun f, double t, double w, double x, double y)
{
   double g = y * t + x;
   return w * y * f(g);
}

double gauss(fun f, double a, double b, int n)
{
   static const struct T
   {
      double abscisass;
      double weights;
   } table[30] =
   {
      /*
       * n = 1
       */
      {               0.0, 2.000000000000000}, /*  0 */
      /*
       * n = 2
       */
      { 0.577350269189626, 1.000000000000000}, /*  1 */
      /*
       * n = 3
       */
      {               0.0, 0.888888888888889}, /*  2 */
      { 0.774596669241483, 0.555555555555556}, /*  3 */
      /*
       * n = 4
       */
      { 0.339981043584856, 0.652145154862546}, /*  4 */
      { 0.861136311594053, 0.347854845137454}, /*  5 */
      /*
       * n = 5
       */
      {               0.0, 0.568888888888889}, /*  6 */
      { 0.538469310105683, 0.478628670499366}, /*  7 */
      { 0.906179845938664, 0.236926885056189}, /*  8 */
      /*
       * n = 6
       */
      { 0.238619186083197, 0.467913934572691}, /*  9 */
      { 0.661209386466265, 0.360761573048139}, /* 10 */
      { 0.932469514203152, 0.171324492379170}, /* 11 */
      /*
       * n = 7
       */
      {               0.0, 0.417959183673469}, /* 12 */
      { 0.405845151377397, 0.381830050505119}, /* 13 */
      { 0.741531185599394, 0.279705391489277}, /* 14 */
      { 0.949107912342759, 0.129484966168870}, /* 15 */
      /*
       * n = 8
       */
      { 0.183434642495650, 0.362683783378362}, /* 16 */
      { 0.525532409916329, 0.313706645877887}, /* 17 */
      { 0.796666477413627, 0.222381034453374}, /* 18 */
      { 0.960289856497536, 0.101228536290376}, /* 19 */
      /*
       * n = 9
       */
      {               0.0, 0.330239355001260}, /* 20 */
      { 0.324253423403809, 0.312347077040003}, /* 21 */
      { 0.613371432700590, 0.260610696402935}, /* 22 */
      { 0.836031107326636, 0.180648160694857}, /* 23 */
      { 0.968160239507626, 0.081274388361574}, /* 24 */
      /*
       * n = 10
       */
      { 0.148874338981631, 0.295524224714753}, /* 25 */
      { 0.433395394129247, 0.269266719309996}, /* 26 */
      { 0.679409568299024, 0.219086362515982}, /* 27 */
      { 0.865063366688985, 0.149451349150581}, /* 28 */
      { 0.973906528517172, 0.066671344308688}, /* 29 */
   };
   double h = 0;
   double x = (b + a) / 2;
   double y = (b - a) / 2;
   switch ( n )
   {
   case 1:
      h += helper(f,  table[0].abscisass, table[0].weights, x, y);
      break;

   case 2:
      h += helper(f,  table[1].abscisass, table[1].weights, x, y);
      h += helper(f, -table[1].abscisass, table[1].weights, x, y);
      break;

   case 3:
      h += helper(f,  table[2].abscisass, table[2].weights, x, y);
      h += helper(f,  table[3].abscisass, table[3].weights, x, y);
      h += helper(f, -table[3].abscisass, table[3].weights, x, y);
      break;

   case 4:
      h += helper(f,  table[4].abscisass, table[4].weights, x, y);
      h += helper(f, -table[4].abscisass, table[4].weights, x, y);
      h += helper(f,  table[5].abscisass, table[5].weights, x, y);
      h += helper(f, -table[5].abscisass, table[5].weights, x, y);
      break;

   case 5:
      h += helper(f,  table[6].abscisass, table[6].weights, x, y);
      h += helper(f,  table[7].abscisass, table[7].weights, x, y);
      h += helper(f, -table[7].abscisass, table[7].weights, x, y);
      h += helper(f,  table[8].abscisass, table[8].weights, x, y);
      break;

   case 6:
      h += helper(f,  table[9].abscisass, table[9].weights, x, y);
      h += helper(f, -table[9].abscisass, table[9].weights, x, y);
      h += helper(f,  table[10].abscisass, table[10].weights, x, y);
      h += helper(f, -table[10].abscisass, table[10].weights, x, y);
      h += helper(f,  table[11].abscisass, table[11].weights, x, y);
      h += helper(f, -table[11].abscisass, table[11].weights, x, y);
      break;

   case 7:
      h += helper(f,  table[12].abscisass, table[12].weights, x, y);
      h += helper(f,  table[13].abscisass, table[13].weights, x, y);
      h += helper(f, -table[13].abscisass, table[13].weights, x, y);
      h += helper(f,  table[14].abscisass, table[14].weights, x, y);
      h += helper(f, -table[14].abscisass, table[14].weights, x, y);
      h += helper(f,  table[15].abscisass, table[15].weights, x, y);
      h += helper(f, -table[15].abscisass, table[15].weights, x, y);
      break;

   case 8:
      h += helper(f,  table[16].abscisass, table[16].weights, x, y);
      h += helper(f, -table[16].abscisass, table[16].weights, x, y);
      h += helper(f,  table[17].abscisass, table[17].weights, x, y);
      h += helper(f, -table[17].abscisass, table[17].weights, x, y);
      h += helper(f,  table[18].abscisass, table[18].weights, x, y);
      h += helper(f, -table[18].abscisass, table[18].weights, x, y);
      h += helper(f,  table[19].abscisass, table[19].weights, x, y);
      h += helper(f, -table[19].abscisass, table[19].weights, x, y);
      break;

   case 9:
      h += helper(f,  table[20].abscisass, table[20].weights, x, y);
      h += helper(f,  table[21].abscisass, table[21].weights, x, y);
      h += helper(f, -table[21].abscisass, table[21].weights, x, y);
      h += helper(f,  table[22].abscisass, table[22].weights, x, y);
      h += helper(f, -table[22].abscisass, table[22].weights, x, y);
      h += helper(f,  table[23].abscisass, table[23].weights, x, y);
      h += helper(f, -table[23].abscisass, table[23].weights, x, y);
      h += helper(f,  table[24].abscisass, table[24].weights, x, y);
      h += helper(f, -table[24].abscisass, table[24].weights, x, y);
      break;

   case 10:
      h += helper(f,  table[25].abscisass, table[25].weights, x, y);
      h += helper(f, -table[25].abscisass, table[25].weights, x, y);
      h += helper(f,  table[26].abscisass, table[26].weights, x, y);
      h += helper(f, -table[26].abscisass, table[26].weights, x, y);
      h += helper(f,  table[27].abscisass, table[27].weights, x, y);
      h += helper(f, -table[27].abscisass, table[27].weights, x, y);
      h += helper(f,  table[28].abscisass, table[28].weights, x, y);
      h += helper(f, -table[28].abscisass, table[28].weights, x, y);
      h += helper(f,  table[29].abscisass, table[29].weights, x, y);
      h += helper(f, -table[29].abscisass, table[29].weights, x, y);
      break;

   default:
      break;
   }
   return h;
}

Again. I have preferred to avoid user input and just generate the output for a number of n cases.

/* my output
n =  1: answer = 78.591673732008658249, error = 2.39883
n =  2: answer = 77.37527840768416354, error = 0.813963
n =  3: answer = 77.064522512218445627, error = 0.409073
n =  4: answer = 76.939340948240840135, error = 0.245971
n =  5: answer = 76.977425880429635185, error = 0.295592
n =  6: answer = 76.840609553794266162, error = 0.117331
n =  7: answer = 76.818123277338983712, error = 0.0880336
n =  8: answer = 76.803122525724845104, error = 0.0684888
n =  9: answer = 76.792617258708546046, error = 0.0548012
n = 10: answer = 76.784974602071045524, error = 0.0448435
*/

jonsca commented: A true, quality post +2

thomas_naveen commented: A good post. Appreciate the time taken +1

Reply to this topic

Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.

jonsca 1,059 Quantitative Phrenologist Team Colleague Featured Poster · Answer 1 · 2010-02-03T22:07:55+00:00

Thanks DS. Can't say I didn't drop the ball on that one. I guess I could say I was going for an overall Gestalt.

sexyzebra19 0 Junior Poster in Training · Answer 2 · 2010-02-03T22:59:36+00:00

Dave - my code runs fine without any errors. It's just that my code actually works for n up to 10 and I didn't want to post it as it's so long, I have obviously forgotten to edit the functions - sorry.

Thanks Jonsca, I appreciate your suggestions, I wasn't sure if repeating code was bad or not, I will try and shorten it down.

sexyzebra19 0 Junior Poster in Training · Answer 3 · 2010-02-03T23:50:12+00:00

OK, posted my full code below - managed to get rid of all those parameters, but I don't see how I can shorten the code, since for each value of N I need to use different abcissass and weight values...is there definitely a way? If there is I will keep thinking :)

#include<iostream>
#include<vector>
#include<cmath>
#include<iomanip>

using namespace std;

typedef double(*fun)(double);

//Declare functions
double gauss(fun f, double &a, double &b, int &n, double &t, double &g, double &h, double &answer, vector<double> &abscisass, vector<double> &weights);
double equation(double t);
double err(double &answer, double &exactanswer, double &error);

//=============================================================================================================

int main()
{
    double a = 0.0;
    double b = 6.0;
    double g = 0.0;
    double h = 0.0;
    double answer = 0.0;
    double exactanswer = 76.750557;
    double error = 0.0;
    double t = 0;
    int n = 0;
    vector<double> abscisass(30,0);
    vector<double> weights(30,0);

    // n = 1
    abscisass[0] = 0.0;
    weights  [0] = 2.000000000000000;
    // n = 2
    abscisass[1] = 0.577350269189626;
    weights  [1] = 1.000000000000000;
    // n = 3
    abscisass[2] = 0.0;
    weights  [2] = 0.888888888888889;
    abscisass[3] = 0.774596669241483;
    weights  [3] = 0.555555555555556;
    // n = 4
    abscisass[4] = 0.339981043584856;
    weights  [4] = 0.652145154862546;
    abscisass[5] = 0.861136311594053;
    weights  [5] = 0.347854845137454;
    // n = 5
    abscisass[6] = 0.0;
    weights  [6] = 0.568888888888889;
    abscisass[7] = 0.538469310105683;
    weights  [7] = 0.478628670499366;
    abscisass[8] = 0.906179845938664;
    weights  [8] = 0.236926885056189;
    // n = 6
    abscisass[9]  = 0.238619186083197;
    weights  [9]  = 0.467913934572691;
    abscisass[10] = 0.661209386466265;
    weights  [10] = 0.360761573048139;
    abscisass[11] = 0.932469514203152;
    weights  [11] = 0.171324492379170;
    // n = 7
    abscisass[12] = 0.0;
    weights  [12] = 0.417959183673469;
    abscisass[13] = 0.405845151377397;
    weights  [13] = 0.381830050505119;
    abscisass[14] = 0.741531185599394;
    weights  [14] = 0.279705391489277;
    abscisass[15] = 0.949107912342759;
    weights  [15] = 0.129484966168870;
    // n = 8
    abscisass[16] = 0.183434642495650;
    weights  [16] = 0.362683783378362;
    abscisass[17] = 0.525532409916329;
    weights  [17] = 0.313706645877887;
    abscisass[18] = 0.796666477413627;
    weights  [18] = 0.222381034453374;
    abscisass[19] = 0.960289856497536;
    weights  [19] = 0.101228536290376;
    // n = 9
    abscisass[20] = 0.0;
    weights  [20] = 0.330239355001260;
    abscisass[21] = 0.324253423403809;
    weights  [21] = 0.312347077040003;
    abscisass[22] = 0.613371432700590;
    weights  [22] = 0.260610696402935;
    abscisass[23] = 0.836031107326636;
    weights  [23] = 0.180648160694857;
    abscisass[24] = 0.968160239507626;
    weights  [24] = 0.081274388361574;
     // n = 10
    abscisass[25] = 0.148874338981631;
    weights  [25] = 0.295524224714753;
    abscisass[26] = 0.433395394129247;
    weights  [26] = 0.269266719309996;
    abscisass[27] = 0.679409568299024;
    weights  [27] = 0.219086362515982;
    abscisass[28] = 0.865063366688985;
    weights  [28] = 0.149451349150581;
    abscisass[29] = 0.973906528517172;
    weights  [29] = 0.066671344308688;

    gauss(equation,a,b,n,t,g,h,answer,abscisass,weights);

    err(answer,exactanswer,error);

    cout << "the answer is " << setprecision(20) << answer << endl;

    cout << "the error is " << setprecision(6) << error << endl;

    return 0;
}

//=============================================================================================================

double gauss(fun f, double &a, double &b, int &n, double &t, double &g, double &h, double &answer, vector<double> &abscisass, vector<double> &weights)
{
    cout << "\n********** GAUSS QUADRATURE **********"<< endl;
    cout << "For the equation f(t) = log(t) + 4t" << endl;
    cout << "Enter the quadrature order, n between 1 and 10 you wish to use: " << endl;
    cin >> n;

    if(n==1)
    {
        t = abscisass[0];
        g = (((b - a)/2)*t) + (b + a)/2;      //coordinate transformation
        h = (weights[0])*((b - a)/2)*(f(g));              // substitution into original equation
        answer = h;
    }
    else if(n==2)
    {
        t = abscisass[1];
        g = (((b - a)/2)*t) + (b + a)/2;      //coordinate transformation
        h = (weights[1])*((b - a)/2)*(f(g));              // substitution into original equation
        t = -1*abscisass[1];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[1])*((b - a)/2)*(f(g));
        answer = h;
    }
    else if(n==3)
    {
        t = abscisass[2];
        g = (((b - a)/2)*t) + (b + a)/2;      //coordinate transformation
        h = (weights[2])*((b - a)/2)*(f(g));              // substitution into original equation
        t = abscisass[3];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[3])*((b - a)/2)*(f(g));
        t = -1*abscisass[3];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[3])*((b - a)/2)*(f(g));
        answer = h;
    }
    else if(n==4)
    {
        t = abscisass[4];
        g = (((b - a)/2)*t) + (b + a)/2;      //coordinate transformation
        h = (weights[4])*((b - a)/2)*(f(g));              // substitution into original equation
        t = -1*abscisass[4];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[4])*((b - a)/2)*(f(g));
        t = abscisass[5];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[5])*((b - a)/2)*(f(g));
        t = -1*abscisass[5];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[5])*((b - a)/2)*(f(g));
        answer = h;
    }
    else if(n==5)
    {
        t = abscisass[6];
        g = (((b - a)/2)*t) + (b + a)/2;      //coordinate transformation
        h = (weights[6])*((b - a)/2)*(f(g));              // substitution into original equation
        t = abscisass[7];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[7])*((b - a)/2)*(f(g));
        t = -1*abscisass[7];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[7])*((b - a)/2)*(f(g));
        t = abscisass[8];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[8])*((b - a)/2)*(f(g));
        t = -1*abscisass[8];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[8])*((b - a)/2)*(f(g));
        answer = h;
    }
    else if(n==6)
    {
        t = abscisass[9];
        g = (((b - a)/2)*t) + (b + a)/2;      //coordinate transformation
        h = (weights[9])*((b - a)/2)*(f(g));              // substitution into original equation
        t = -1*abscisass[9];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[9])*((b - a)/2)*(f(g));
        t = abscisass[10];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[10])*((b - a)/2)*(f(g));
        t = -1*abscisass[10];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[10])*((b - a)/2)*(f(g));
        t = abscisass[11];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[11])*((b - a)/2)*(f(g));
        t = -1*abscisass[11];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[11])*((b - a)/2)*(f(g));
        answer = h;
    }
    else if(n==7)
    {
        t = abscisass[12];
        g = (((b - a)/2)*t) + (b + a)/2;      //coordinate transformation
        h = (weights[12])*((b - a)/2)*(f(g));              // substitution into original equation
        t = abscisass[13];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[13])*((b - a)/2)*(f(g));
        t = -1*abscisass[13];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[13])*((b - a)/2)*(f(g));
        t = abscisass[14];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[14])*((b - a)/2)*(f(g));
        t = -1*abscisass[14];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[14])*((b - a)/2)*(f(g));
        t = abscisass[15];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[15])*((b - a)/2)*(f(g));
        t = -1*abscisass[15];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[15])*((b - a)/2)*(f(g));
        answer = h;
    }
    else if(n==8)
    {
        t = abscisass[16];
        g = (((b - a)/2)*t) + (b + a)/2;      //coordinate transformation
        h = (weights[16])*((b - a)/2)*(f(g));              // substitution into original equation
        t = -1*abscisass[16];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[16])*((b - a)/2)*(f(g));
        t = abscisass[17];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[17])*((b - a)/2)*(f(g));
        t = -1*abscisass[17];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[17])*((b - a)/2)*(f(g));
        t = abscisass[18];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[18])*((b - a)/2)*(f(g));
        t = -1*abscisass[18];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[18])*((b - a)/2)*(f(g));
        t = abscisass[19];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[19])*((b - a)/2)*(f(g));
        t = -1*abscisass[19];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[19])*((b - a)/2)*(f(g));
        answer = h;
    }
        else if(n==9)
    {
        t = abscisass[20];
        g = (((b - a)/2)*t) + (b + a)/2;      //coordinate transformation
        h = (weights[20])*((b - a)/2)*(f(g));              // substitution into original equation
        t = abscisass[21];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[21])*((b - a)/2)*(f(g));
        t = -1*abscisass[21];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[21])*((b - a)/2)*(f(g));
        t = abscisass[22];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[22])*((b - a)/2)*(f(g));
        t = -1*abscisass[22];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[22])*((b - a)/2)*(f(g));
        t = abscisass[23];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[23])*((b - a)/2)*(f(g));
        t = -1*abscisass[23];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[23])*((b - a)/2)*(f(g));
        t = abscisass[24];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[24])*((b - a)/2)*(f(g));
        t = -1*abscisass[24];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[24])*((b - a)/2)*(f(g));
        answer = h;
    }
    else if(n==10)
    {
        t = abscisass[25];
        g = (((b - a)/2)*t) + (b + a)/2;      //coordinate transformation
        h = (weights[25])*((b - a)/2)*(f(g));              // substitution into original equation
        t = -1*abscisass[25];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[25])*((b - a)/2)*(f(g));
        t = abscisass[26];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[26])*((b - a)/2)*(f(g));
        t = -1*abscisass[26];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[26])*((b - a)/2)*(f(g));
        t = abscisass[27];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[27])*((b - a)/2)*(f(g));
        t = -1*abscisass[27];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[27])*((b - a)/2)*(f(g));
        t = abscisass[28];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[28])*((b - a)/2)*(f(g));
        t = -1*abscisass[28];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[28])*((b - a)/2)*(f(g));
        t = abscisass[29];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[29])*((b - a)/2)*(f(g));
        t = -1*abscisass[29];
        g = (((b - a)/2)*t) + (b + a)/2;
        h = h + (weights[29])*((b - a)/2)*(f(g));
        answer = h;
    }
}

//=============================================================================================================

double equation (double t)
{
    return log(t) +4*t;
}

double err(double &answer, double &exactanswer, double &error)
{
    error = fabs(((answer - exactanswer)/exactanswer)*100);
}

sexyzebra19 0 Junior Poster in Training · Answer 4 · 2010-02-05T03:15:54+00:00

Dave, I've just had a quick read through your post...you are actually amazing. Thank you SO much for taking the time to do that, just from glancing over it I feel that I've learned a lot, your advice is brilliant and you've explained it so well. I'm now off to study your post in detail, I cannot thank you enough :)

sexyzebra19 0 Junior Poster in Training · Answer 5 · 2010-02-05T06:21:23+00:00

Running into another problem now which I've been looking at for hours and can't seem to fix.

I'm trying to improve the method implemented, by instead of evaluating over integral 0 and 6, I introduce the variable M, which splits the integral into M sections.

So for the minute I'm using M = 2, which should evaluate the function at 0 and 3, and then 3 and 6, which I'm hoping will be more accurate.

I tried do do this below, but I'm getting NaN for all answers and really can't see why.

The problem seems to be here:

double h = 0;
   double h1 = (b-a)/m;
   double x = (b2 + a2) / 2;
   double y = (b2 - a2) / 2;

as I left it as

double h = 0;
   double h1 = (b-a)/m;
   double x = (b + a) / 2;
   double y = (b - a) / 2;

and it ran OK but gave me the wrong answer.

Can you Dave, or anybody else see what is causing this?

#include <iostream>
#include <cmath>
#include <iomanip>
using namespace std;

typedef double(*fun)(double);

double equation(double t)
{
   return log(t) + 4*t;
}

double err(double answer, double exactanswer){
   return fabs(((answer - exactanswer) / exactanswer) * 100);
}

double pattern(fun f, double t, double w, double x, double y)
{
   double g = y * t + x;
   return w * y * f(g);
}

double gauss(fun f, double a, double b, int n, int m);

int main()
{
   double exactanswer = 76.750557;
   for ( int n = 1; n <= 10; ++n )
   {
      double answer = gauss(equation, 0, 6, n, 2);
      cout << "n = " << setw(2) << n
           << ": answer = " << setprecision(20) << answer
           << ", error = "  << setprecision(6)  << err(answer, exactanswer)
           << endl;
   }
   return 0;
}



double gauss(fun f, double a, double b, int n, int m)
{
   static const struct T
   {
      double abscisass;
      double weights;
   }
   table[30] =
   {
      // n = 1
      {               0.0, 2.000000000000000}, // 0
      // n = 2
      { 0.577350269189626, 1.000000000000000}, // 1
      // n = 3
      {               0.0, 0.888888888888889}, // 2
      { 0.774596669241483, 0.555555555555556}, // 3
      // n = 4
      { 0.339981043584856, 0.652145154862546}, // 4
      { 0.861136311594053, 0.347854845137454}, // 5
      // n = 5
      {               0.0, 0.568888888888889}, // 6
      { 0.538469310105683, 0.478628670499366}, // 7
      { 0.906179845938664, 0.236926885056189}, // 8
      // n = 6
      { 0.238619186083197, 0.467913934572691}, // 9
      { 0.661209386466265, 0.360761573048139}, // 10
      { 0.932469514203152, 0.171324492379170}, // 11
      // n = 7
      {               0.0, 0.417959183673469}, // 12
      { 0.405845151377397, 0.381830050505119}, // 13
      { 0.741531185599394, 0.279705391489277}, // 14
      { 0.949107912342759, 0.129484966168870}, // 15
      // n = 8
      { 0.183434642495650, 0.362683783378362}, // 16
      { 0.525532409916329, 0.313706645877887}, // 17
      { 0.796666477413627, 0.222381034453374}, // 18
      { 0.960289856497536, 0.101228536290376}, // 19
      // n = 9
      {               0.0, 0.330239355001260}, // 20
      { 0.324253423403809, 0.312347077040003}, // 21
      { 0.613371432700590, 0.260610696402935}, // 22
      { 0.836031107326636, 0.180648160694857}, // 23
      { 0.968160239507626, 0.081274388361574}, // 24
      // n = 10
      { 0.148874338981631, 0.295524224714753}, // 25
      { 0.433395394129247, 0.269266719309996}, // 26
      { 0.679409568299024, 0.219086362515982}, // 27
      { 0.865063366688985, 0.149451349150581}, // 28
      { 0.973906528517172, 0.066671344308688}, // 29
   };

   double h = 0;
   double h1 = (b-a)/m;
   double x = (b2 + a2) / 2;
   double y = (b2 - a2) / 2;

   switch ( n )
   {
   case 1:
   for (int j = 1; j <= m; j++)
        {
            a2 = a + (j-1)*h1;
            b2 = a + j*h1;
            h += pattern(f,  table[0].abscisass, table[0].weights, x, y);
        }
      break;

   case 2:
      h += pattern(f,  table[1].abscisass, table[1].weights, x, y);
      h += pattern(f, -table[1].abscisass, table[1].weights, x, y);
      break;

   case 3:
      h += pattern(f,  table[2].abscisass, table[2].weights, x, y);
      h += pattern(f,  table[3].abscisass, table[3].weights, x, y);
      h += pattern(f, -table[3].abscisass, table[3].weights, x, y);
      break;

   case 4:
      h += pattern(f,  table[4].abscisass, table[4].weights, x, y);
      h += pattern(f, -table[4].abscisass, table[4].weights, x, y);
      h += pattern(f,  table[5].abscisass, table[5].weights, x, y);
      h += pattern(f, -table[5].abscisass, table[5].weights, x, y);
      break;

   case 5:
      h += pattern(f,  table[6].abscisass, table[6].weights, x, y);
      h += pattern(f,  table[7].abscisass, table[7].weights, x, y);
      h += pattern(f, -table[7].abscisass, table[7].weights, x, y);
      h += pattern(f,  table[8].abscisass, table[8].weights, x, y);
      break;

   case 6:
      h += pattern(f,  table[9].abscisass, table[9].weights, x, y);
      h += pattern(f, -table[9].abscisass, table[9].weights, x, y);
      h += pattern(f,  table[10].abscisass, table[10].weights, x, y);
      h += pattern(f, -table[10].abscisass, table[10].weights, x, y);
      h += pattern(f,  table[11].abscisass, table[11].weights, x, y);
      h += pattern(f, -table[11].abscisass, table[11].weights, x, y);
      break;

   case 7:
      h += pattern(f,  table[12].abscisass, table[12].weights, x, y);
      h += pattern(f,  table[13].abscisass, table[13].weights, x, y);
      h += pattern(f, -table[13].abscisass, table[13].weights, x, y);
      h += pattern(f,  table[14].abscisass, table[14].weights, x, y);
      h += pattern(f, -table[14].abscisass, table[14].weights, x, y);
      h += pattern(f,  table[15].abscisass, table[15].weights, x, y);
      h += pattern(f, -table[15].abscisass, table[15].weights, x, y);
      break;

   case 8:
      h += pattern(f,  table[16].abscisass, table[16].weights, x, y);
      h += pattern(f, -table[16].abscisass, table[16].weights, x, y);
      h += pattern(f,  table[17].abscisass, table[17].weights, x, y);
      h += pattern(f, -table[17].abscisass, table[17].weights, x, y);
      h += pattern(f,  table[18].abscisass, table[18].weights, x, y);
      h += pattern(f, -table[18].abscisass, table[18].weights, x, y);
      h += pattern(f,  table[19].abscisass, table[19].weights, x, y);
      h += pattern(f, -table[19].abscisass, table[19].weights, x, y);
      break;

   case 9:
      h += pattern(f,  table[20].abscisass, table[20].weights, x, y);
      h += pattern(f,  table[21].abscisass, table[21].weights, x, y);
      h += pattern(f, -table[21].abscisass, table[21].weights, x, y);
      h += pattern(f,  table[22].abscisass, table[22].weights, x, y);
      h += pattern(f, -table[22].abscisass, table[22].weights, x, y);
      h += pattern(f,  table[23].abscisass, table[23].weights, x, y);
      h += pattern(f, -table[23].abscisass, table[23].weights, x, y);
      h += pattern(f,  table[24].abscisass, table[24].weights, x, y);
      h += pattern(f, -table[24].abscisass, table[24].weights, x, y);
      break;

   case 10:
      h += pattern(f,  table[25].abscisass, table[25].weights, x, y);
      h += pattern(f, -table[25].abscisass, table[25].weights, x, y);
      h += pattern(f,  table[26].abscisass, table[26].weights, x, y);
      h += pattern(f, -table[26].abscisass, table[26].weights, x, y);
      h += pattern(f,  table[27].abscisass, table[27].weights, x, y);
      h += pattern(f, -table[27].abscisass, table[27].weights, x, y);
      h += pattern(f,  table[28].abscisass, table[28].weights, x, y);
      h += pattern(f, -table[28].abscisass, table[28].weights, x, y);
      h += pattern(f,  table[29].abscisass, table[29].weights, x, y);
      h += pattern(f, -table[29].abscisass, table[29].weights, x, y);
      break;

   default:
      break;
   }
   return h;
}

jonsca 1,059 Quantitative Phrenologist Team Colleague Featured Poster · Answer 6 · 2010-02-05T06:36:13+00:00

Where are a2 and b2 supposed to come from? You attempt to write to them below those equations in case 1 but they are never declared anywhere (not even in your global scope -- which is good that you didn't, but then they'd be somewhere!).

sexyzebra19 0 Junior Poster in Training · Answer 7 · 2010-02-05T07:04:31+00:00

Where are a2 and b2 supposed to come from? You attempt to write to them below those equations in case 1 but they are never declared anywhere (not even in your global scope -- which is good that you didn't, but then they'd be somewhere!).

Sorry sorry, I did actually declare them within the gauss function originally but I was cutting and pasting them into different sections trying to get it to work and obviously forgot to put them in :-/

I had:

double a2 = 0;
   double b2= 0;
   double h = 0;
   double h1 = (b-a)/m;
   double x = (b2 + a2) / 2;
   double y = (b2 - a2) / 2;

when I tested it, but it still seems to find a problem with that...

jonsca 1,059 Quantitative Phrenologist Team Colleague Featured Poster · Answer 8 · 2010-02-05T10:13:18+00:00

You only have 4 calls to your helper function in your case 5, btw.

From what it looks like, you're basically setting x = 0, (since (0 + 0)/2), y = 0 (from (0-0)/2 )and those values are never changing as you pass them in with the different abscissas. So why not just call all those methods with x = 0,y= 0. I'm not sure the way that you have it is exactly what you want.