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

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

Yehh, I finally made the program in C to calculate roots of a user entered polynomial in a user entered interval, it gives output - the multiplicity of the root proceeded by the root . The algorithm used is that firstly we divide the interval into degree+1 number of equally spaced intervals then Newton-Raphson method is applied in each interval to calculate maximum number of roots. Several checks are also included for the cases when user gives invalid input.For some polynomials it might give undesired output

#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("\tYou entered an alphabet. Kindly, retry with valid entry\n");
    else
        printf("\tYou entered a character. Kindly, retry with valid entry\n");

    return 0;
}

char signum(float x)
{
    if (x>=0) return '+';
    else return '-';
}



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

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

  return fc;
}


int fact(int n)
{
  if (n>0) return n*fact(n-1);
  else return 1;
}



double D(int k,int a,double x)               //calculates the derivative of the function required in the formmula in main
{
  if(x!=0 && a>=k)  return fact(a)*pow(x,a-k)/fact(a-k);
  else return 0;  
}


double derv(double *q,int degree,double x,int order)               //calculates the derivative of the function required in the formmula in main
{
  double dv=0.0;
  int i;
  for(i=0;i<=degree;i++)
  {
    dv=dv+(*q)*D(order,degree-i,x);
    q++;
  }
  return dv;
}



int main()
{
    double coeff[1000],t,x1=0.0,roots[100];
    int degree=-1,count=0,i,j=0,k,l=0,flag=1;
    char ch;
    float a,b;
    coeff[0]=0;
    printf("\t**This Programme helps you find roots of a polynomial using a manipulated version of Newton-Raphson Method in a real interval [a,b]**\n");


    while (degree<0)
    {
        printf("\tEnter a NON-NEGATIVE INTEGRAL value for the DEGREE of the polynomial: ");
        while(scanf("%d",&degree)!=1)
        {
            printf("\tInvalid Entry. Kindly, retry with valid entry: ");
            while(getchar()!='\n'){}
        }
}

int w[degree];
printf("\tEnter the coefficients...\n");

while(coeff[0]==0 && degree!=0)
{
    printf("\t-->for x^%d (Remember that this coefficient cannot be zero since degree is %d): ", degree,degree);
    while(scanf("%lf",&coeff[0])!=1)
    {
        check(coeff[0]);
    }
}

if (degree==0)
{
    printf("\t-->for x^%d: ",0);
    while(scanf("%lf",&coeff[0])!=1)
    {
        check(coeff[0]);
    }
}

for(i=1;i<=degree;i++)
{
    printf("\t-->for x^%d: ",degree-i);
    while(scanf("%lf",&coeff[i])!=1)
    {
        check(coeff[i]);
    }
}

  while(getchar()!='\n'){}                                 //this clears the buffer accumulated upto pressing enter
    printf("\tThe input polynomial is:\n");
  printf("\t");
  for(i=0;i<=degree;i++)
  {
    printf("%c %.2fx^%d ",signum(coeff[i]),fabs(coeff[i]),degree-i);
  }

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

  /* fflush(stdout);
  fflush(stdin);*/
  printf("\tFor the interval domain [a,b]; enter 'a' followed by 'b': ");
  scanf("%f %f",&a,&b);


  while(getchar()!='\n'){}
    printf("\tSo, you entered [%f,%f].\n",a,b);
  //while(getchar()!='\n'){}
  fflush(stdout);
  fflush(stdin);


  double d=(b-a)/(degree+1), val[degree+2],x=a;
  for (k = 0; k<=degree+1; k++)
  {
    val[k]= a+d*k;
  //    printf("%lf\t", val[k]);
  }

  while(l<(degree+2))
  {

    if(derv(&coeff[0],degree,x,1)<=0.00001)
    {
        if(func(&coeff[0],degree,x)<=0.00001)
        {
            w[j]=2;
            roots[j]=x;
        //printf("Axk before %f\n", x);
            j++;
        //  printf("jC is %d\n", j);
            l++;
            x=val[l];   // 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
         //printf("Axk after %f\n", x);
            t=0.0;
            x1=0.0;
            count=(-1);
        }
        else
        {
            t=0.0;
            x1=0.0;
            count=(-1);
            l++;
          //printf("Bxk before%f\n", x);
            x=val[l];
          //  printf("Bxk after%f\n", x);
        }
    }
       else
       {
          t=x1;
          x1=(x-(func(&coeff[0],degree,x)/derv(&coeff[0],degree,x,1)));           //this is the newton's formula
       // printf("jB is %d\n", j);
          //printf("Cxk before%f\n", x);
          x=x1;
          //printf("f %f\td %f\tCxk after%f\tc %i\n",func(&coeff[0],degree,x), derv(&coeff[0],degree,x,1), x,count);
          if(count>500)
          {

            if(j>degree)
                break;
              printf("\tPolynomial started to diverge. So, no roots can be found\n");
            l++;
            //printf("Dxk before %f\n", x);
            x=val[l];
            //printf("Dxk after %f\n", x);
            t=0.0;
            x1=0.0;
            count=(-1);

          }     
          if((fabs(t - x1))<=0.00001 && count!=0)
          {
            w[j]=1;
            if(j>degree)
                break;  

            roots[j]=x;

            j++;
            l++;
        //  printf("jC is %d\n", j);
            x=val[l];   // 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
         //printf("Exk after%f\n", x);
            t=0.0;
            x1=0.0;
            count=(-1);
            flag=0;
            if(derv(&coeff[0],degree,x,1)<=0.00001) w[j]=2;
            else w[j]=1;
           }
           count++;
       } 
   }

   if(flag==0)
   {
    printf("\tThe roots are = \t");
     // int p=0;
    for(j=0;j<degree;j++)
    {

        if(j==0 && roots[0]>=a && roots[0]<=b)
        {
            printf(" %d %f\n", w[0], roots[0]);
       //p++;
        }

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


    }
   } 
   else
    printf("\tNo roots found in the interval you entered.\n");

   return 0;

}

Answer 2

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.