VDC117
VDC117

Reputation: 3

Logical error in bisection algorithm code

I am writing a program to demonstrate the bisection algorithm(from numerical methods).

What I am doing is this:

  1. defined a function F(int), which takes the integer and returns the value of the polynomial at that integer.

  2. In the bisection() function:

    a) a,b are the initial approximations.

    b) the for() construct finds the value of a after which the sign(magnitude) of F(a) changes.

    c) the printfs are for troubleshooting purposes(except for the last one).

    d) int prec is for the precision required for the solution (no. of iterations), where no. of iterations, N=log((b-a)/E), E=10^(-prec).

What I'm getting is this:

Output:

a=1     F(a)=-0.207936

a=1, b=2
Approximation 0: x = 0  a=1072693248   b=1
Approximation 1: x = 0  a=1072693248   b=1
Approximation 2: x = 0  a=1072693248   b=1
Approximation 3: x = 0  a=1072693248   b=1
Approximation 4: x = 0  a=1072693248   b=1

The solution is: x = 1.000000

I have tried commenting the N=... statement and assigning constant integer to N, but to no effect.

Interestingly, commenting out all statements with variable x in the for() construct does not alter values of a and b.

Oddly, if the statement x=(a+b)/2; is commented out, then the value of a is affected by the value with which x is initialized:

i) a=1072693248, when x=1,

ii)when x=0, a=0,b=0

iii)when x=-1, a=-1074790400, b=1

I am using Microsoft Visual C++ 2010 Express for compilation.

Please help me out. Here is the code:

#include<stdio.h>
#include<math.h>

#define MIN -10

double F(int x)
{
    double y=(x*x)-(2*x)+0.792064;
//  printf("\nx=%d y=%d",x,y);
    return y;
}

void bisection(int prec)
{
    int a,b=0;
    int c,N=10,i;
    double x=0;

    for(a=MIN;(F(a)*F(a-1))>0;a++);
    printf("\na=%d  F(a)=%f",a,F(a));


    b=a+1;
    printf("\n\na=%d, b=%d",a,b);
    N=(log((float)(b-a))+5);

    for(i=0;i<N;i++)
    {
            x=(a+b)/2;
            printf("\nApproximation %d: x = %d\ta=%d   b=%d",i,x,a,b);
            if((F(a)*F(x)>0))
                a=x;
            else
                b=x;
    }
    printf("\n\nThe solution is: x = %f",x);
    getchar();
}

int main()
{   
    bisection(4);
    return 0;
}

Upvotes: 0

Views: 161

Answers (1)

Sufian Latif
Sufian Latif

Reputation: 13356

It's because a and b are declared as ints. When you write x = (a + b) / 2;, all the elements on the right side are integers, so they will be evaluated to an integer. Moreover, as the difference of a and b is 1, the value would be equal to the smaller one of them, which is a. That's why the value of x becomes 1.0. To fix it, the upper and lower limits should be declared as doubles, not ints. When a and b are integers, their values cannot converge to the actual value of x.

Also, you're using %d in printf to write a double value. It should be %lf.

Upvotes: 2

Related Questions