[Newbie] What's wrong with this code?

  • Thread starter Manuele Santoprete
  • Start date
M

Manuele Santoprete

I wrote this code basically the trapezoidal rule to compute integrals.
If I compile it with different versions of gcc and run the executable I
get different results. In some cases I get 0.335 that seems the right
result, in other cases I get something close to 2.5.

Any help would be appreciated.

Thanks!!


#include <stdio.h>
#include <math.h>
# define N 11
double Q[]={0,0.01,0.04,0.09,0.16,0.25,0.36,0.49,0.64,0.81,1},
dt[]={0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1};
/*double L[N], h[N];*/
double integral;
double trapzd (double L[N], double h[N-2]);
double S;


int
main(void)
{
integral=trapzd(Q,dt);
printf("%g\n",integral);
return 0;
}

double trapzd (double L[],double h[])
{

int i;
for(i=0;i<=N-2;i++)
{
S=S+0.5*h*(L[i+1]+L);
}
return S;
}
 
X

xhungab

/* -------------------------------------
Hello,

Can I offer my code,
I am not able to correct your code.

Thank
------------------------------------- */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* -------------------------------------
------------------------------------- */
double f(
double x)
{
return(sqrt(x));
}
char feq[] = "sqrt(x)";
/* -------------------------------------
------------------------------------- */
double F(
double x)
{
return( (2./3.) * pow(x,(3./2.)));
}
char Feq[] = " (2./3.)* x**(3./2.)";
/* -------------------------------------
------------------------------------- */
double trapezoid(
double (*P_f)(double x),
double a,
double b,
int n
)
{
int i = 0;
double m = 0.;
double M = 0.;

for(i = 0; i <= n; i++)
{
if(i ==0 || i== n){m = 1.;}
else {m = 2.;}

M += m * (*P_f)(a + i*(b-a)/n);
}

return( ((b -a)*M) / (2*n) );
}
/* ------------------------------- MAIN */
int main(void)
{
double (*P_f)(double x);

int n = 2*30000;
double a = 0.;
double b = 1.;
double M = 0.;
/*----------------------------- PROGRAM */
P_f = f;

printf(" With the trapezoidal's rule.\n\n");
M = trapezoid((*P_f),a,b,n);
printf(" (%.3f\n", b);
printf(" int( (%s) * dx = %.12f\n", feq, M);
printf(" (%.3f\n\n\n\n", a);

printf(" With antiderivative of f.\n\n");
printf(" (%.3f\n", b);
printf(" int( (%s) * dx = %.12f\n", feq, F(b)-F(a));
printf(" (%.3f\n\n\n\n", a);

printf("\n Press return to continue");
fflush(stdout);
getchar();

return 0;
}
 
A

Ancient_Hacker

Manuele Santoprete wrote:

I don't see anywhere that you initialize "S". Maybe you want to set it
to zero at the beginning of the function? And maybe also move it
inside the function.
 
M

msantopr

Tanks to everybody!
I think Ancient_Hacker is right I didn't initialize S! I should have
set it to zero as he suggested.

Thanks!
 
D

dcorbit

Tanks to everybody!
I think Ancient_Hacker is right I didn't initialize S! I should have
set it to zero as he suggested.

Thanks!

That's nonsense. S is global and hence is initialized to zero.

The problem lies elsewhere.

When you fix something that is not broken and the problem goes away, it
isn't fixed.
It looks like a compiler bug to me, but I have not examined it very
carefully.
 
D

dcorbit

Manuele said:
I wrote this code basically the trapezoidal rule to compute integrals.
If I compile it with different versions of gcc and run the executable I
get different results. In some cases I get 0.335 that seems the right
result, in other cases I get something close to 2.5.

Any help would be appreciated.

Thanks!!


#include <stdio.h>
#include <math.h>
# define N 11
double Q[]={0,0.01,0.04,0.09,0.16,0.25,0.36,0.49,0.64,0.81,1},
dt[]={0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1};
/*double L[N], h[N];*/
double integral;
double trapzd (double L[N], double h[N-2]);
double S;


int
main(void)
{
integral=trapzd(Q,dt);
printf("%g\n",integral);
return 0;
}

double trapzd (double L[],double h[])
{

int i;
for(i=0;i<=N-2;i++)
{
S=S+0.5*h*(L[i+1]+L);
}
return S;
}


/*
I would rework it to something more like this:
*/

#include <stdlib.h> /* for size_t */
double trapzd(const double y[], double h, size_t length)
{
double result = 0;
size_t i;

for (i = 0; i < length - 1; i++) {
result += 0.5 * h * (y[i + 1] + y);
}
return result;
}

#ifdef UNIT_TEST
#include <stdio.h>
#include <float.h>

static const double y[] = {0, 0.01, 0.04, 0.09, 0.16, 0.25,
0.36, 0.49, 0.64, 0.81, 1};

int main(void)
{
double integral = trapzd(y, 0.1, sizeof y / sizeof y[0]);
printf("%*.*g\n", DBL_DIG+4, DBL_DIG, integral);
return 0;
}
#endif

/*
I guess that the problems you encountered came from using global
variables.

If you use global variables and do not reset them between runs, then
the values will accumulate and you will get wrong results.

C:\tmp>cl /W4 /Ox /DUNIT_TEST trap.c
Microsoft (R) 32-bit C/C++ Optimizing Compiler Version 14.00.50727.42
for 80x86
Copyright (C) Microsoft Corporation. All rights reserved.

trap.c
Microsoft (R) Incremental Linker Version 8.00.50727.42
Copyright (C) Microsoft Corporation. All rights reserved.

/out:trap.exe
trap.obj

C:\tmp>trap
0.335
*/
 
A

Ancient_Hacker

That's nonsense. S is global and hence is initialized to zero.

The problem lies elsewhere.

When you fix something that is not broken and the problem goes away, it
isn't fixed.
It looks like a compiler bug to me, but I have not examined it very
carefully.

It's usually a poor idea to write a function that can only be called
once.
 
D

dcorbit

Ancient_Hacker wrote:
[snip]
It's usually a poor idea to write a function that can only be called
once.

I usually call main() only once in a program.
;-)

Really, even if I intend to call a function only once, I will break
things down into small blocks that are easy to understand and make
these little jobs into functions (which is a side point and not really
related to your concept -- which I tend to agree with).
If I only intend to call it once (and within my program), I will
probably make it static.

In this case, the problem was in the use of global variables, none of
which needed to be global. Even the test array could easily have been
auto, but that one did not hurt anything.
 
B

Barry Schwarz

I wrote this code basically the trapezoidal rule to compute integrals.
If I compile it with different versions of gcc and run the executable I
get different results. In some cases I get 0.335 that seems the right
result, in other cases I get something close to 2.5.

My system says .335 also and Excel agrees. You could strategically
place a few printf statements to see where the erroneous value comes
from.
Any help would be appreciated.

Thanks!!


#include <stdio.h>
#include <math.h>
# define N 11
double Q[]={0,0.01,0.04,0.09,0.16,0.25,0.36,0.49,0.64,0.81,1},
dt[]={0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1};
/*double L[N], h[N];*/
double integral;
double trapzd (double L[N], double h[N-2]);
double S;


int
main(void)
{
integral=trapzd(Q,dt);
printf("%g\n",integral);
return 0;
}

double trapzd (double L[],double h[])
{

int i;
for(i=0;i<=N-2;i++)
{
S=S+0.5*h*(L[i+1]+L);
}
return S;
}



Remove del for email
 

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments. After that, you can post your question and our members will help you out.

Ask a Question

Members online

No members online now.

Forum statistics

Threads
473,744
Messages
2,569,484
Members
44,903
Latest member
orderPeak8CBDGummies

Latest Threads

Top