Horner's Method

Horner's method is an efficient way of evaluating polynomials and their derivatives at a given point. It is also used for a compact presentation of the long division of a polynomial by a linear polynomial.

A polynomial P_n(x) can be written in the descending order of the powers of x:

Pn(x) = anxn + an-1xn-1 + ... + a1x + a0

or in the ascending order of the exponents:

Pn(x) = a0 + a1x + ... + an-1xn-1 + anxn.

The Horner scheme computes the value P_n(z) of the polynomial P at x = z as the sequence of steps starting with the leading coefficient a_n using:

	bk = z·bk+1 + ak
	ck = z·ck+1 + bk

                The next value of the x_n in the iterative step is found by,

                x_k+1 = x_k – (b₀/c₁), where b₀ = P_n(x_k) = b₀ and c₁ = P’_n(x_k) in shortcut.

Generally, the iteration is continued until the required amount of precision is obtained.

Example: Solve the equation x³+9x²-18 using Horner’s Method where z= x₀ = 1.5.

            Here, a₀=-18, a₁=0, a₂=9, and a₃=1.

            Starting from a₃, b₃=c₃=a₃=1.

            Hence,

            b₂ = 9+1x1.5 = 10.5                          c₂ = 10.5+1x1.5 = 12

            b₁ = 0+10.5x1.5 = 15.75                  c₁ = 15.75+12x1.5 = 33.75

            b₀ = -18+15.75x1.5 = 5.625                 

            So, x₁ = x₀ – (b₀/c₁) = 1.5 – (5.625/33.75) = 1.333

            Þ x₁ = z = 1.333

            

            Now, for the new value of z,

            b₂ = 9 + 1x1.333 = 10.333                               c₁ = 28.881

            b₁ = 0+10.333x1.333 = 13.774

            b₀ = -18 + 13.774x1.333 = 0.361

            So, x₂ = x₁ – (b₀/c₁) = 1.333 – (0.361/28.881) = 1.320

            Þ x₂ = z = 1.320

            Again, for new value of z,

            b₀ = -0.018           c₁ = 28.987

            So, x₃ = x₂ – (b₀/c₁) = 1.320 – (-0.018/28.987) = 1.320

Therefore, Approximate Root = 1.320 correct to 3 decimal places. 

C program for Cubic Spline Interpolation

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

void main(){
clrscr();
char choice='y';
int n,i,j,k;
float h[10],a,b,c,d,sum,s[10]={0},x[10],F[10],f[10],p,m[10][10]={0},temp;

printf("No of samples? ");
scanf("%d",&n);
printf("\nEnter all sample points: ");
for(i=0;i<n;i++)
scanf("%f%f",&x[i],&f[i]);
for(i=n-1;i>0;i--){
F[i]=(f[i]-f[i-1])/(x[i]-x[i-1]);
h[i-1]=x[i]-x[i-1];
}
//*********** formation of h, s , f matrix **************//
for(i=1;i<n-1;i++){
m[i][i]=2*(h[i-1]+h[i]);
if(i!=1){
m[i][i-1]=h[i-1];
m[i-1][i]=h[i-1];
}
m[i][n-1]=6*(F[i+1]-F[i]);
}
//*********** forward elimination **************//
for(i=1;i<n-2;i++){
temp=(m[i+1][i]/m[i][i]);
for(j=1;j<=n-1;j++)
m[i+1][j]-=temp*m[i][j];
}
//*********** backward substitution *********//
for(i=n-2;i>0;i--){
sum=0;
for(j=i;j<=n-2;j++)
sum+=m[i][j]*s[j];
s[i]=(m[i][n-1]-sum)/m[i][i];
}
while(choice=='y'){
printf("\nEnter x : ");
scanf("%f",&p);
for(i=0;i<n-1;i++)
if(x[i]<=p&&p<=x[i+1]){
a=(s[i+1]-s[i])/(6*h[i]);
b=s[i]/2;
c=(f[i+1]-f[i])/h[i]-(2*h[i]*s[i]+s[i+1]*h[i])/6;
d=f[i];
sum=a*pow((p-x[i]),3)+b*pow((p-x[i]),2)+c*(p-x[i])+d;
}
printf("\nCoefficients of sub intervals are : %f\n%f\n%f\n%f",a,b,c,d);
printf("\nFunctional value is: %f",sum);
printf("\nContinue (y/n) ? ");
scanf("%c",&choice);
}
getch();
}

Sample Output:

No of samples? 5

Enter all sample points: 1     2

2     3

3     4

4     5

5     6

Enter x : 3.5

Coefficients of sub intervals are : 0.000000

0.000000

1.000000

4.000000

Functional value is: 4.500000

Continue (y/n) ?

©Dixit Bhatta 2013

C++ program for Newton-Raphson Method

#include<iostream.h>
#include<conio.h>
#include<math.h>
#define EPS 0.000001
#define MAXIT 20
#define F(x) (x)*(x)*(x)+(x)*(x)-3*(x)-3
#define FD(x) 3*(x)*(x)+2*(x)-3

void main(){
clrscr();
int count;
float x0, xn, fx, fdx, fxn;
cout<<"             NEWTON RAPHSON'S METHOD";
   cout<<endl<<"x^3+x^2-3x-3";
cout<<endl<<"Enter the initial value of x:";
cin>>x0;
count=1;

begin:
fx=F(x0);
fdx=FD(x0);
xn=x0-fx/fdx;
if (fabs((xn-x0)/xn) < EPS){
cout<<"Approximate Root: "<<xn<<endl;
fxn =F(xn);
cout<<"Functional Value: "<<fxn<<endl;
cout<<"No. of Iterations: "<<count<<endl;
}
else{
x0=xn;
count=count+1;
if(count<MAXIT){
goto begin;
}
else{
cout<<"SOLUTION DOES NOT CONVERGE!!"<<endl;
cout<<"Iterations: "<<MAXIT<<endl;
}
}
getch();
}

Sample Output:

                    NEWTON RAPHSON’S METHOD

x^3+x^2-3x-3

Enter the initial value of x:2

Approximate Root: 1.73205

Functional Value: -2.94213e-07

No. of Iterations: 4

©Dixit Bhatta 2013

C Program for Bisection Method

Basic Theory:

The bisection method is a method to find the real root of an equation in which the interval is always divided into half for every new iterative step. If a function changes sign over an interval, the functional value at the mid-point is evaluated. The location of the root is then determined to be within the sub-interval where the sign change occurs. The sub-interval thus becomes the interval for next iteration. The process is repeated until the root is known to the required precision.

Following is a numerical example of Bisection method.

Consider finding the root of f(x) = x² - 3. Let Error = 0.01 and start with the interval [1, 2].

Table 1: Bisection method applied to f(x) = x² - 3.

a	b	f(a)	f(b)	c = (a + b)/2	f(c)	Exchange	Error
1.0	2.0	-2.0	1.0	1.5	-0.75	a = c	0.5
1.5	2.0	-0.75	1.0	1.75	0.062	b = c	0.25
1.5	1.75	-0.75	0.0625	1.625	-0.359	a = c	0.125
1.625	1.75	-0.3594	0.0625	1.6875	-0.1523	a = c	0.0625
1.6875	1.75	-0.1523	0.0625	1.7188	-0.0457	a = c	0.0313
1.7188	1.75	-0.0457	0.0625	1.7344	0.0081	b = c	0.0156
1.7188	1.7344	-0.0457	0.0081	1.7266	-0.0189	a = c	0.0078

Thus, with the seventh iteration, we note that the final interval, [1.7266, 1.7344], has a width less than the value of Error and |f(1.7344)| < 0.01, and therefore we chose b = 1.7344 to be our approximation of the root.

Program:
#include<stdio.h>

#include<conio.h>

#include<math.h>

#define f(x) x*x*x-2*x-5 //define the function here

void main(){

clrscr();

float x1 = 2, x2 =3, E=0.0005, x0; //define initial values and errors here

float f1 = f(x1), f2 = f(x2), f0;

float root;

                if(f1*f2>0){

                                goto last;

                }

printf("\n\t\t\tBISECTION METHOD");

printf("\n\n  FINDING ROOT OF x^3-2x-5");

printf("\n\n  X1\t  X2\t  X0\tf(X0)\t  Error");

printf("\n  -------------------------------------");

back:

x0 = (x1+x2)/2;

f0 = f(x0);

printf("\n  %.3f\t%.3f\t%.3f\t%.3f\t%.5f",x1,x2,x0,f(x0),fabs((x2-x1)/x2));

//print in tabular form

                if(f1*f0<0){

                                x2= x0;

                }

                else{

                                x1=x0;

                                f1=f0;

                }

                if(fabs((x2-x1)/x2)<E){

                                root = (x1+x2)/2; //display root if required precision is met

                                printf("\n\n  Approximate Root = %.5f", root);

                                goto last;

                }

                else{

                                goto back;

                }

last:

getch();

}

Sample Output:

     BISECTION METHOD
FINDING ROOT OF x^3-2x-5

X1          X2       X0      f(X0)      Error
------------------------------------------------
2.000   3.000   2.500   5.625    0.33333
2.000   2.500   2.250   1.891    0.20000
2.000   2.250   2.125   0.346    0.11111
2.000   2.125   2.062   -0.351   0.05882
2.062   2.125   2.094   -0.009   0.02941
2.094   2.125   2.109   0.167    0.01471
2.094   2.109   2.102   0.079    0.00741
2.094   2.102   2.098   0.035    0.00372
2.094   2.098   2.096   0.013    0.00186
2.094   2.096   2.095   0.002    0.00093
2.094   2.095   2.095   0.002    0.00047


Approximate Root = 2.09424 

©Dixit Bhatta 2013

C Program for Lagrange's Interpolation

#include<stdio.h>

#include<conio.h>

#define MAX 10

void main(){

clrscr();

float x[MAX+1], fx[MAX+1], num, den, val, y=0;

int i, j, n;

printf("\t\tLAGRANGE'S INTERPOLATION");

printf("\n\nEnter the degree of interpolation(n): ");

scanf("%d", &n);

for(i=0; i<=n; i++){

printf("\n x%d: ",i);

scanf("%f", &x[i]);

printf("\n fx%d: ",i);

scanf("%f", &fx[i]);

}

printf("\nEnter x to estimate fx: ");

scanf("%f", &val);

for (i=0; i<=n; i++){

num=1;

den=1;

for (j=0; j<=n; j++)

if(j!=i){

num *= val-x[j];

den *= x[i]-x[j];

}

y+=(num/den)*fx[i];

}

printf("\nApproximate value at %.3f: %.3f",val,y);

getch();

}

Sample Output:

LAGRANGE’S INTERPOLATION
Enter the degree of interpolation(n): 2
x0: 1
fx0: 1
x1: 2
fx1: 1.4142
x2: 3
fx2: 1.7321

Enter x to estimate fx: 2.5
Approximate value at 2.500: 1.585

©Dixit Bhatta 2013

C Program for Horner's Method

#include<stdio.h>

#include<conio.h>

#include<math.h>

float p[6], ply[6], q[6];

float synth(int m, float r){

                int i;

                q[0] = p[0];

                for(i=1;i<=m;i++){

                                q[i] = (q[i-1]*r)+p[i];

                }

                return(q[m]);

}

void main(){

                clrscr();

                int m,i,flag=0;

                float r, x,x1, b0, c1;

                printf("\t\tHORNER'S METHOD");

                printf("\n Enter the highest degree of the equation (max 5): ");

                scanf("%d",&m);

                                for(i=0;i<=m;i++){

                                printf("\n Coefficient x[%d] = ",m-i);

                                scanf("%f",&p[i]);

                                ply[i] = p[i];

                                }

                printf("\n Enter the initial value x0 : ");

                scanf("%f",&r);

                x = r;

                do{

                   printf("\n%f\n",x);

                   b0 = synth(m,x);

                                for(i=0;i<=5;i++){

                                                p[i]=q[i];

                                }

                   c1 = synth(m-1,x);

                   printf("\n");

                   printf("\tb0 = %f\tc1= %f",b0,c1);

                   x1 = x - (b0/c1);

                   if(fabs(x1-x) <= 0.0009){

                                flag = 1;

                   }

                   x = x1;

                                for(i=0;i<=5;i++){

                                                p[i]=ply[i];

                                }

                   getch();

                }while(flag!=1);

                printf("\nApproximate root = %f", x1);

getch();

}

Sample Output:

                          HORNER’S METHOD
Enter the highest degree of the equation (max 5): 3
Coefficient X[3] : 1
Coefficient X[2] : 9
Coefficient X[1] : 0
Coefficient X[0] : -18
Enter the initial value X0 : 1.5
1.500000
              b0 = 5.625000   c1= 33.750000
1.333333
              b0 = 0.370371   c1= 29.333332
1.328787
              b0 =  0.002071   c1= 29.005529
Approximate root = 1.320636

©Dixit Bhatta 2013

C Program for Birge-Vieta Method

Basic Theory:

The Birge-Vieta method investigates the real root of a polynomial, say x_r, by assuming that (x-x_r) is a factor of that polynomial. Starting with an initial value or approximation x₀ to x_r, we use some iterative method to improve the value of the approximation such that R(x_r)®0. It makes the use of the Newton-Raphson method for iteration in order to improve the value of x_r.

x_n+1=x_n – [f(x_n)/f ’(x_n)]

We deploy the Synthetic Division method in order to find out R(x_r) in this method. The remainder of first division gives the functional value while the remainder second division gives the value of derivative. We continue the iteration until the result of desired precision is obtained.

Example:  Solve the equation x³+x²-3x-3 using Birge-Vieta Method where x₀ = 2.

Using the synthetic division,

2|1         1              -3            -3

  |           2              6              6

  |1         3              3              3¬f(x₀)

  |           2              10

  |1         5              13¬f ’(x₀)

Now, x₁ = 2 – 3/13 = 1.7692

1.7692|1              1              -3            -3

           |                 1.7692   4.8992   3.3602

           |1              2.7692   1.8992   0.3602¬f(x₁)

           |                 1.7692   8.0292

           |1              4.5384   9.9285¬f ’(x₁)

Now, x₂ = 1.7692 – 0.3602/9.9285 = 1.7329

1.7329|1              1              -3            -3

           |                 1.7329   4.7358   3.0079

           |1              2.7329   1.7358   0.0079¬f(x₂)

           |                 1.7329   7.7387

           |1              4.4658   9.4745¬f ’(x₂)

Now, x₃ = 1.7329 – 0.0079/9.4745 = 1.7320

Therefore, Approximate Root = 1.732 correct to 3 decimal places.

Program:

#include<stdio.h>

#include<conio.h>

#include<math.h>

float p[6], ply[6],q[6];

float synth(int m, float r){

                int i;

                q[0] = p[0];

                for(i=1;i<=m;i++){

                                q[i] = (q[i-1]*r)+p[i];

                }

                printf("\n");

                for(i=0;i<m;i++){

                                printf("\t%f",q[i]);

                }

                printf("\t%f",q[m]);

                return(q[m]);

}

void main(){

                clrscr();

                int m,i,flag=0;

                float r, x,x1, fx, fdx;

                printf("\t\tBIRGE-VIETA METHOD");

                printf("\n Enter the highest degree of the equation (max 5): ");

                scanf("%d",&m);

                                for(i=0;i<=m;i++){

                                printf("\n Coefficient x[%d] = ",m-i);

                                scanf("%f",&p[i]);

                                ply[i] = p[i];

                                }

                printf("\n Enter the initial value x0 : ");

                scanf("%f",&r);

                x = r;

                do{

                   printf("\n%f\n",x);

                   fx = synth(m,x);

                                for(i=0;i<=m;i++){

                                                p[i]=q[i];

                                }

                   fdx = synth(m-1,x);

                   x1 = x - (fx/fdx);

                                if(fabs(x1-x) <= 0.0009){

                                                flag = 1;

                                 }

                   x = x1;

                                for(i=0;i<=5;i++){

                                                p[i]=ply[i];

                                }

                }while(flag!=1);

                printf("\nApproximate root = %f", x1);

getch();

}

Sample Output:

                        BIRGE-VIETA METHOD
Enter the highest degree of the equation (max 5): 3
Coefficient x[3] = 1
Coefficient x[2] = 1
Coefficient x[1] = -3
Coefficient x[0] = -3
Enter the initial value x0 : 2
2.000000
              1.000000   3.000000   3.000000   3.000000
              1.000000   5.000000   13.000000
1.769231
              1.000000   2.769231   1.899408   0.360492
              1.000000   4.538462   9.928994
1.732924
              1.000000   2.732924   1.735948   0.008266
              1.000000   4.465847   9.474921
Approximate root = 1.732051

©Dixit Bhatta 2013

C Program for Synthetic Division with algorithm and example

Algorithm of Synthetic Division:

Given a polynomial of form p(x) = a_nxⁿ + a_n-1x^n-1+…+ a₁x+ a₀, we can divide it by a linear factor x-r, where ‘r’ is a constant, using following steps.

a.       Write the r on the left side of a vertical bar and the coefficients in decreasing degree.

r| a_n               a_n-1         …….        a₁            a₀

b.      Pass the coefficient of the highest degree below the horizontal line as it is

 r| a_n              a_n-1         …….        a₁            a₀

  |.                                                                   .

  | a_n             

c.       Multiply the passed coefficient with ‘r’ and add to the next  coefficient

 r| a_n              a_n-1         …….        a₁            a₀

  |.                 a_n                                                  .

  | a_n          a_n-1+ a_n

d.      Continue the step ‘c’ until all coefficients are covered.

e.      The numbers below a_n ……. a₁ give the coefficients of the quotient whose degree is 1 less than the polynomial, the number below a₀ is the remainder.

Example of Synthetic Division:  x² + 5x+ 6 by x – 1

Writing the coefficients of the polynomial and dividing by the linear factor.

1|1         5              6

  |           1              6

  |1         6              12

Therefore, the quotient Q(x) = x + 6 and remainder R(x) = 12.

Program:

#include<stdio.h>

#include<conio.h>

void main(){

                clrscr();

                int poly[6], m, r, i, q[6];

                printf("\t\tSYNTHETIC DIVISION");

                printf("\n Enter the highest degree of the equation (max 5): ");

                scanf("%d",&m);

                for(i=0;i<=m;i++){

                                printf("\n Coefficient x[%d] = ", m-i);

                                scanf("%d",&poly[i]);

                }

                printf("\n Enter the value of constant in (x-r) : ");

                scanf("%d",&r);

                q[0] = poly[0];

                for(i=1;i<=m;i++){

                                q[i] = (q[i-1]*r)+poly[i];

                }

                printf("\n Coefficients of quotient are: \n");

                for(i=0;i<m;i++){

                                printf("\t%d",q[i]);

                }

                printf("\n Remainder is: %d", q[m]);

getch();

}

Sample Output:

SYNTHETIC DIVISION
Enter the highest degree of the equation (max 5): 2

Coefficient x[0] = 1

Coefficient x[1] = 5

Coefficient x[2] = 6

Enter the value of constant in (x—r) : 1

Coefficients of quotient are:
1        6
Remainder is: 12

©Dixit Bhatta 2013

C Program for Secant Method

#include<stdio.h>

#include<conio.h>

#include<math.h>

#define E 0.0001

#define F(x) x*x - 4*x - 10

void main(){

  float x1,x2,x3,f1,f2,t;

  clrscr();

  printf("\nEnter the value of x1: ");

  scanf("%f",&x1);

  printf("\nEnter the value of x2: ");

  scanf("%f",&x2);

  printf("\n______________________________________________\n");

  printf("\n    x1\t  x2\t  x3\t     f(x1)\t f(x2)");

  printf("\n______________________________________________\n");

  do{

                f1=F(x1);

                f2=F(x2);

                x3=x2-((f2*(x2-x1))/(f2-f1));

                printf("\n%f   %f   %f   %f   %f",x1,x2,x3,f1,f2);

                x1=x2;

                x2=x3;

                                if(f2<0)

                                                t=fabs(f2);

                                else

                                                t=f2;

  }while(t > E);

printf("\n______________________________________________\n");

printf("\n\nApp.root = %f",x3);

getch();

}

Sample Output:

Enter the value of x1: 4
Enter the value of x2: 2
-----------------------------------------------
x0            x1        x2           f(x0)           f(x1)
-----------------------------------------------
4.000   2.000   9.0000   -10.000    -14.000
2.000   9.000   4.0000   -14.000    35.000
9.000   4.000   5.1111   35.000     -10.000
4.000   5.111   5.9565   -10.000    -4.321
5.111   5.957   5.7225   -4.321      1.654
5.957   5.722   5.7411   1.654       -0.143
5.722   5.741   5.7417   -0.143     -0.004
5.741   5.742   5.7417   -0.004      0.000
-----------------------------------------------

App.root = 5.741657

©Dixit Bhatta 2013