C++ simpson's method -loop doesn't stop

[ggeo1]

c++ simpson's method --loop doesn't stop

Hello, i have done the simpson's method and it works fine.My problem is that it doesn't stop when it gives me the desired result but the loop continues for ever..

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <iomanip>
#include <cmath>

using namespace std;

double f(double x){
double y;
y=sqrt(x-0.1);
return y;

}

double simpson (double down,double up,int n){

    double h,sumeven,sumodd,all,result;

    h=(up-down)/n;
     sumeven=0.0;
    for (int i=1;i<=n/2.0-1.0;i++){

        sumeven+=f(down+2.0*i*h);//coefficients

    }
    sumodd=0.0;
    for (int i=1;i<=n/2.0;i++){
     sumodd+=f(down+h*(2.0*i-1.0));//coefficients

    }


    all=4*sumodd+2*sumeven;

    result=(h/3.0)*(f(down)+f(up)+all);
    return result;

}


int main()
{
    double eps=1e-6;//accuracy
    double exact=1.1767695;//exact solution for the integral
    double error=1.0;
    double result;

    int n=1;//initial point
    result=simpson(1.0,2.0,n);

    while (fabs(error-exact)>eps) {
        result=simpson(1.0,2.0,n);
        cout <<"\nFor n = "<<n<<",error = "<<fabs(error-exact)<<",value = "<<result;

    n++;
    }

    return 0;
}

I have done other methods too like that and they stop,but this one not.

Thanks!

 

[jtbell]

In your main function, you define

    double eps=1e-6;//accuracy
    double exact=1.1767695;//exact solution for the integral
    double error=1.0;

but then in your while-loop you don't change any of those variables, so the condition 'fabs(error-exact)>eps' is always true.

 

[Mark44]

This loop looks flaky to me.

for (int i=1;i<=n/2.0-1.0;i++){
        sumeven+=f(down+2.0*i*h);//coefficients
}

In the first iteration of the while loop in main, n is 1, so the test expression in the for loop above is false, so the for loop doesn't execute.

When n == 1, n/2.0 - 1.0 is 0.5 - 1.0 == -.5, and i <= -.5 is false. You should not have to compare the counter in the for loop (i) with a double.

 

[ggeo1]

Hello,you are right but

if i do :

 int n=1;//initial point
    result=simpson(1.0,2.0,n);

    while (fabs(result-exact)>eps) {
        result=simpson(1.0,2.0,n);
        cout <<"\nFor n = "<<n<<",error = "<<fabs(result-exact)<<",value = "<<result;

    n++;
    }

then the loop stops at step n=8 ,but the problem is that at step n=2 it gives me the right solution and then goes over again.

 

[ggeo1]

Ok, i changed the for loops : for (int i=1;i<=n/2-1;i++)
for (int i=1;i<=n/2;i++)
(it was a mistake ,i didn't want to insert double values in the loops!)

But now again , it gives me the exit as a say above..

 

[Mark44]

I think that you are not using Simpson's Rule correctly. Each time you call simpson(), the 3rd argument should be an even integer.

In each iteration of Simpson's Rule the interval [a, b] should be broken up into an even number of subintervals.

 

[ggeo1]

I used the same code in mathematica and it works fine.It gives me the right answer after 2 iterations.The same is happening here but it continues to step 8.(it gives me the right result every "even" steps)..

 

[Mark44]

Instead of incrementing n by 1 in each while loop iteration, try doubling n each time.

So instead of doing this: n++;
Do this: n *= 2;

See if that makes a difference.

You've changed your code since the first post. Can you show us what you have now?
Also, it would help to see the exact output.

 

[ggeo1]

Here is my code :

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <iomanip>
#include <cmath>

using namespace std;

double f(double x){
double y;
y=sqrt(x-0.1);
return y;

}

double simpson (double down,double up,int n){

    double h,sumeven,sumodd,all,result;

    h=(up-down)/n;
     sumeven=0.0;
    for (int i=1;i<=n/2-1;i++){

        sumeven+=f(down+2.0*i*h);//coefficients

    }
    sumodd=0.0;
    for (int i=1;i<=n/2;i++){
     sumodd+=f(down+h*(2.0*i-1.0));//coefficients

    }


    all=4*sumodd+2*sumeven;

    result=(h/3.0)*(f(down)+f(up)+all);
    return result;

}


int main()
{
    double eps=1e-6;//accuracy
    double exact=1.1767695;//exact solution for the integral
    double error=1.0;
    double result;

    int n=1;//initial point
   

    while (fabs(result-exact)>eps) {
        result=simpson(1.0,2.0,n);
        cout <<"\nFor n = "<<n<<",error = "<<fabs(result-exact)<<",value = "<<result;

    n++;
    }

    return 0;
}

and the exit is :

For n = 1,error = 0.401073,value = 0.775696
For n = 2,error = 0.000110833,value = 1.17666
For n = 3,error = 0.424624,value = 0.752146
For n = 4,error = 8.23034e-06,value = 1.17676
For n = 5,error = 0.263326,value = 0.913444
For n = 6,error = 1.67268e-06,value = 1.17677
For n = 7,error = 0.190651,value = 0.986119
For n = 8,error = 5.17684e-07,value = 1.17677

If i use n*=2 instead of n++ ,it gives :

For n = 1,error = 0.401073,value = 0.775696
For n = 2,error = 0.000110833,value = 1.17666
For n = 4,error = 8.23034e-06,value = 1.17676
For n = 8,error = 5.17684e-07,value = 1.17677

which is the same but with step 2 of course.