Like Tree1Likes
  • 1 Post By dmh2000

Double integral help

This is a discussion on Double integral help within the C Programming forums, part of the General Programming Boards category; I think I'm doing everything correctly but for some reason the program refuses to compute, it just freezes. I'm given ...

  1. #1
    Registered User
    Join Date
    Aug 2011
    Posts
    8

    Double integral help

    I think I'm doing everything correctly but for some reason the program refuses to compute, it just freezes.

    I'm given an equation for a curve, and for certain values of x and y along that I want to take the double integral and set that over the double integral of the whole curve and compute that fraction.

    To do this I loop through values of x and y and then multiply them to get the double integral of the whole curve.

    For the isolated section of the curve I substitute a new variable z for y and then within the x loop I cycle through z values for each value of x and if x and z are within certain values I compute the smaller integral.

    I then set the two integrals over each other and compute that fraction.

    What am I doing wrong?

    Code:
    #include <stdio.h>
    #include <math.h>
    #define STEP 5.556e-5
    
    int main() {
    float a=0.3489, b=0.255, R=0.68085, xlimit=0.2799, zlimit=-0.0008; 
    float utotalarea=0; 
    float vtotalarea=0; 
    float x, xarea, xtotalarea=0; 
    float y, yarea, ytotalarea=0; 
    float z, zarea; 
    float result, result2, fraction;
    //First loop through x values.
    for (x=STEP; x<=10; x=x+STEP){
    	xarea=0.5*(exp(-log(2)*pow(x/a,2)-1)+exp(-log(2)*pow((x-STEP)/a,2)-1))*STEP;
    	xtotalarea=xtotalarea+xarea;
    		//Next, loop through z values for each x value. Note: z is really y, its a different variable here in order to calculate the integral for certain x and y values.  
    		for (z=STEP; z<=10; z=z+STEP){
    			zarea=0.5*(exp(-log(2)*pow(z/b,2)-1) + exp(-log(2)*pow((z-STEP)/b,2)-1))*STEP;
    			//If x and z are within certain values take integral.			
    			if((pow(x-xlimit,2)+pow(z-zlimit,2))<=pow(R,2)){
    				utotalarea=utotalarea+xarea;				
    				vtotalarea=vtotalarea+zarea;
    			}
    		}
    	}
    //Now, loop through all y values.
    for (y=STEP; y<=10; y=y+STEP){
    	yarea=0.5*(exp(-log(2)*pow(y/b,2)-1)+exp(-log(2)*pow((y-STEP)/b,2)-1))*STEP;
    	ytotalarea=ytotalarea+yarea;
    	}
    //Compute double integral for the whole curve.
    result=xtotalarea*ytotalarea;
    //Compute double integral for the certain x and z values.
    result2=utotalarea*vtotalarea;
    //Set the smaller integral over the larger integral.
    fraction=result2/result;
    //print the fraction. 
    printf("%.5f\n", fraction);
    }

  2. #2
    - - - - - - - - oogabooga's Avatar
    Join Date
    Jan 2008
    Posts
    2,808
    Since both your outer and inner loop are repeating 10 / .00005556 = 180000 times, your inner loop is repeating about 32.4 billion times. It might finish in a two or three hours (days?).

  3. #3
    Registered User
    Join Date
    Mar 2011
    Posts
    546
    for kicks just increase STEP to 1 and see if it runs. the answer won't be what you want but you will see if it runs to completion.
    Salem likes this.

  4. #4
    Registered User
    Join Date
    Sep 2008
    Location
    Toronto, Canada
    Posts
    1,832
    I would suggest to never use float. Use double instead. float is a relic from decades ago where one would sacrifice numerical accuracy for memory. Nowadays all math.h functions use and return double. It's actually faster to stick with double than implicit casting to float all the time.

  5. #5
    Algorithm Dissector iMalc's Avatar
    Join Date
    Dec 2005
    Location
    New Zealand
    Posts
    6,308
    Quote Originally Posted by nonoob View Post
    I would suggest to never use float. Use double instead. float is a relic from decades ago where one would sacrifice numerical accuracy for memory. Nowadays all math.h functions use and return double. It's actually faster to stick with double than implicit casting to float all the time.
    Only when memory bandwidth is not a concern. Afterall, 16-bit halffloats are popular in 3D graphics because they take half the memory of floats, which is great if you have a lot of them, say for a z-buffer.
    Last time I checked the difference in speed on floats vs doubles in my software 3D engine, floats were still faster.

    Having said that, I agree that doubles are best in this case.
    My homepage
    Advice: Take only as directed - If symptoms persist, please see your debugger

    Linus Torvalds: "But it clearly is the only right way. The fact that everybody else does it some other way only means that they are wrong"

  6. #6
    Registered User
    Join Date
    Jun 2005
    Posts
    6,529
    I'd suggest it's the algorithm (the complete logic of how you organise the calculations) that needs refining. Changing from float to double might paper over the gaps, but can equally fail again if some parameters get tweaked.

    However, if you're not willing to apply some thought into changing algorithm, you can simplify the code you have......

    It is generally considered a no-no to use a floating point variable to control a loop. Floating point variables are NOT real numbers, they are a finite approximation, so all sorts of errors creep in. Particularly if you are trying to iterate (i.e. add a fixed value each time through the loop). Thanks to finite precision of floating point, rounding error can build up, the value of the control variable will be affected by error, and may well be too inaccurate for needs. Do something like.
    Code:
       for (int i = 1; i <= some_number; ++i)
       {
              x = STEP * i;
               .....
    rather than using x as the control variable and adding a small value of STEP to it.

    You might want to look into simplifying various expressions in the loops, to avoid replication of calculations.

    x*x is better - computationally - for squaring x, than pow(x,2).

    The fact you have exp(-log(2) .... in several expressions suggests the formula you're using actually seeks to raise 2 to various powers. In the calculation of xarea, for example, exp(-log(2)*pow(x/a,2)-1) can be simplified as 1.0/(pow(2.0, x*x/a*a)*exp(1.0)). That can be simplified even further, since a and exp(1.0) are both non-varying values in your code. So you might hoist a fair few repetitive calculations out of the loops (only compute them once, rather than for every nested iteration, and store their values in a temporary variable).

    As I said, however, you would be better off looking at whether there is a better algorithm. There are lots of libraries out there for doing numerical integrations of various special functions.
    Last edited by grumpy; 12-22-2011 at 03:27 PM.
    Right 98% of the time, and don't care about the other 3%.

    If I seem grumpy in reply to you, it is likely you deserve it. Suck it up, sunshine, and read this, this, and this before posting again.

  7. #7
    Captain Crash brewbuck's Avatar
    Join Date
    Mar 2007
    Location
    Portland, OR
    Posts
    7,263
    I agree that for this problem, using float doesn't make a lot of sense. But the memory saving aspects of float are certainly not irrelevant even with huge main memories these days. You might have 10 GB of RAM in your box, but the poor little core processing your data still only has 32k of L1 cache (64k on AMD)
    Code:
    //try
    //{
    	if (a) do { f( b); } while(1);
    	else   do { f(!b); } while(1);
    //}

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Calculating Integral; Please help!
    By sam... in forum C Programming
    Replies: 4
    Last Post: 02-23-2011, 09:20 AM
  2. Calculating an integral
    By kamme in forum C Programming
    Replies: 5
    Last Post: 12-04-2010, 04:00 PM
  3. integral counterpart to pow()?
    By cyberfish in forum C++ Programming
    Replies: 5
    Last Post: 11-23-2007, 10:21 PM
  4. integral promotions
    By St0rM-MaN in forum C Programming
    Replies: 7
    Last Post: 06-28-2007, 05:37 PM
  5. integral of 1/x
    By Silvercord in forum A Brief History of Cprogramming.com
    Replies: 43
    Last Post: 03-19-2004, 07:47 AM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21