Newton's method program (in C) loop running infinitely

Question

This code(attached with the post) in C uses Newton - Raphson method to find roots of a polynomial in a particular interval.

This code works perfectly fine for some polynomials like x^3 + x^2 + x + 1 but runtime gets infinite for some polynomials like x^3 - 6*x^2 + 11*x - 6 . That is this code works fine for polynomials having one or zero root in the entered interval but if more than one roots are present then it runs indefinitely.

Please let me know if someone finds a solution for this. I have written comments in the code to guide the reader but still if anyone finds it difficult to understand, can ask in the comments, I'll explain it.

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

int check(float num)                          //just a function to check for  the correct input
{
    char c;
    scanf("%c",&c);

    if(isalpha((int)c))
        printf("you entered an alphabet\n");

    else
        printf("you entered a character, please retry\n");

    return 0;
}

float func(float *p,int order,double x)                     //calculates the value of the function required in the formmula in main
{
    double fc=0.0;
    int i;

    for(i=0;i<=order;i++)
    {
        fc=fc+(double)(*p)*pow(x,(double)i);
        p++;
    }

    return fc;
}

float derv(float *q,int order,double x)               //calculates the derivative of the function required in the formmula in main
{     
    double dv=0.0,i;

    for(i=1;i<=order;i++)
    {
        dv=dv+(double)(*q)*(pow(x,(double)(i-1)))*(double)i;
        q++;
    }

    return dv;
}


int main()
{
    float coeff[1000];
    int order,count=0,i,j=0;
    char ch;
    float a,b;
    double val[5];

    printf("roots of polynomial using newton and bisection method\n");
    printf("enter the order of the equation\n");

    while(scanf("%d",&order)!=1)
    {
        printf("invalid input.please retry\n");
        while(getchar()!='\n'){}          
    }     

    printf("enter the cofficients\n");

    for(i=0;i<=order;i++)
    {
        printf("for x^%d  :",order-i);
        printf("\n");

        while(scanf("%f",&coeff[i])!=1)
        {
            check(coeff[i]);
        }   
    }

    while(getchar()!='\n'){}                                 //this clears the buffer accumulated upto pressing enter

    printf("the polynomial you entered is :\n");

    for(i=0;i<=order;i++)
    {
        printf(" %fx^%d ",coeff[i],order-i);
    }

    printf("\n");

    //while(getchar()!='\n'){};

    /* fflush(stdout);
    fflush(stdin);*/

    printf("plese enter the interval domain [a,b]\n");
    printf("enter a and b:\n");
    scanf("%f %f",&a,&b);

    while(getchar()!='\n'){}

    printf("the entered interval is [%f,%f]",a,b);

    fflush(stdout);
    fflush(stdin);

    //this array is used to choose a different value of x to apply newton's formula recurcively in an interval to scan it roperly for 3 roots

    val[0]=(double)b;       
    val[1]=(double)a;
    val[2]=(double)((a+b)/2);

    double t,x=val[0],x1=0.0,roots[10];

    while(1)
    {

        t=x1;
        x1=(x-(func(&coeff[0],order,x)/derv(&coeff[0],order,x)));           //this is the newton's formula

        x=x1;

        if((fabs(t - x1))<=0.0001 && count!=0)
        {
            roots[j]=x;
            j++;
            x=val[j];   // every time a root is encountered this stores the root in roots array and runs the loop again with different value of x to find other roots
            t=0.0;
            x1=0.0;
            count=(-1);

            if(j==3)
                break;
        }

        count++;
    }

    printf("the roots are = \n");

    int p=0;

    for(j=0;j<3;j++)
    {
        if(j==0 && roots[j]>=a && roots[j]<=b)
        {
            printf("  %f  ",roots[j]);
            p++;
        }

        if(fabs(roots[j]-roots[j-1])>0.5 && j!=0 && roots[j]>=a && roots[j]<=b)
        {
            printf(" %f  ",roots[j]);
            p++;
        }
    }

    if(p==0)
        printf("Sorry,no roots are there in this interval \n");

    return 0;
}

Answer 1

You're not calculating the function or the derivative correctly, because you store the coefficients in reverse order but you aren't accounting for that.

When you print out the equation, you do account for it, by printing order-i:

printf(" %fx^%d ",coeff[i],order-i);

So you need to do the same thing in func:

fc=fc+(double)(*p)*pow(x,(double)(order-i));

and derv:

dv=dv+(double)(*q)*(pow(x,(double)((order-i)-1)))*(double)(order-i);

The reason it's working for polynomials like x^3 + x^2 + x + 1 is because in that example all the coefficients are the same so it won't make a difference if you read the array forwards or backwards.

In addition, you may need to account for other reasons the method can fail to converge, as mentioned by Johnathon Leffler in a comment. You can set a maximum number of iterations for your loop and just break out if it exceeds the maximum.

A good way to debug something like this (besides using a debugger of course) is to add a few extra printf statements to show the values that are being calculated. You can sanity check the output by entering the equation into a Google search, it will give an interactive graph of the function.