NUMERICAL METHODS
Programs
1. Successive approximation method
Algorithm
1.
2.
3.
4.
5.
6.
Start
Define the iterative equation
Input x0, initial guess
Input the accuracy
Input the max iterations suggested, n
For i =1 to n, increment in steps of 1
a.
b.
c.
d.
x1= f(x0)
Error=|(x1-x0)/x1|
x0 =x1
If (error<=accuracy) output root=x1 and go to step 8
7. Output - roots does not converge in n iterations
8. Stop
1. Successive approximation method
void main()
{ float f(float );
float x0,x1, error, acc;
int i, n;
cin>>x0>>acc>>n;
for(i=0; i<n; i++)
{
x1=f(x0);
error = fabs((x1-x0)/x1);
x0 = x1;
if (error <= acc)
{
cout<< x1;
go to out;
}
}
cout<<root does not converge;
out:
}
float f(float x)
{
float y=pow(x, 5) - 2 * x * x + 3 * x - 1;
return(y);
}
2. Newton Raphson method
1. Start
2. Input x0, acc
3. K=0
4. x[k+1] = x[k] - f(x[k]) / f1(x[k])
5. Error=fabs ((x[k+1] - x[k]) / x[k+1])
6. If (error>= accy)
k=k+1; go to step 4;
7. Print x[k+1]
8. stop
2. Newton Raphson method
void main()
{
double x[200], accy;
int k = 0;
double f(double a);
double f1(double a);
void where(void);
where();
cout <<"\nEnter the initial guess, x0.\n";
cin >>x[0];
cout <<"Enter the accuracy required.\n";
cin >>accy;
x[1] = x[0] - f(x[0]) / f1(x[0]);
2. Newton Raphson method
cout <<k <<"\t" <<x[0] <<"\t<<f(x*0+) <<"\t<<f1(x*0+) <<"\t" <<x[1]
<<"\t<<fabs((x[1] - x[0]) / x[1]) <<"\n";
while(fabs((x[k+1] - x[k]) / x[k+1]) >= accy)
{
k++;
x[k+1] = x[k] - f(x[k]) / f1(x[k]);
cout <<k <<"\t" <<x[k] <<"\t"
<<f(x[k]) <<"\t" <<f1(x[k])
<<"\t<<x*k+1+ <<"\t<<fabs((x[k+1] - x[k]) / x[k+1]) <<"\n";
}
cout <<"One of the real roots is: " <<x[k+1];
2. Newton Raphson method
// Subroutine to evaluate the given function
double f(double a)
{
double y;
y = a * a * a - 3 * a * a + a + 1;
return(y);
}
// Subroutine to evaluate the derivative of the given function
double f1(double a)
{
double y;
y = 3 * a * a - 6 * a + 1;
return(y);
}
2. Newton Raphson method
//Subroutine to determine the intervals where the roots lie
void where(void)
{
float y0, y1, j;
for(j = -50; j <= 50; j++)
{
y0 = f(j);
y1 = f(j + 1);
if((y0 * y1) < 0)
cout <<"A root exists between " <<j <<" and " <<j + 1 <<". \n";
else
continue;
}
return;
}
3. Regula Falsi method
1.
Read x0, x1, error, n
2.
f0 =f(x0)
3.
f1= f(x1)
4.
For (i = 1 to n) in steps of 1 do
5.
x2= (x0*f1-x1*f0)/(f1-f0)
6.
f2= f(x2)
7.
If |f2|<= error, then write convergent soln, exit
8.
If (f0 *f2 >0) x0=x2
else
End for loop
9. stop
x1=x2
3. Regula Falsi method
void main()
{ int n, i;
float x0,x1,x2,error;
cin>>x0>>x1;
cin>>error;
float f0, f1, f2;
f0=f(x0);
f1=f(x1);
do
{
x2=(x0*f1- x1*f0)/(f1- f0);
f2=f(x2);
if ((f0*f2)>0)
{
x0 = x2;
f0 = f2;
}
else
{
x1=x2;
f1 = f2;
}
} while (f2>error);
}
4. Fitting a straight line
void main()
{
float x[20], y[20];
int i, n;
float s1, s2, s3, s4;
cout <<"\nEnter the no. of
values.\n";
cin >>n;
for( i = 0; i < n; i++)
{
cout <<"Enter x[" <<i+1 <<"]
and y[" <<i+1 <<"]";
cin >>x[i] >>y[i];
}
s1 = s2 = s3 = s4 = 0;
for(i = 0; i < n; i++)
{
s1 += x[i];
s2 += x[i] * x[i];
s3 += y[i];
s4 += x[i] * y[i];
}
float d = n * s2 - s1 * s1;
float a = (s2 * s3 - s1 * s4) / d;
float b = (n * s4 - s1 * s3) / d;
cout <<"The equation for the
straight line is: " <<a <<"x + "
<<b;
}
5. Fitting a parabola
void main()
{
float x[20], y[20], a[5][5];
int i, n;
float s1, s2, s3, s4, s5, s6, s7;
void guass(float a[5][5]);
cout <<"\nEnter the no. of values.\n";
cin >>n;
for( i = 0; i < n; i++)
{
cout <<"Enter x[" <<i+1 <<"]
and y[" <<i+1 <<"]";
cin >>x[i] >>y[i];
}
s1 = s2 = s3 = s4 = s5 = s6 = s7 = 0;
for(i = 0; i < n; i++)
{
s1 += x[i];
s2 += y[i];
s3 += x[i] * x[i];
s4 += x[i] * x[i] * x[i];
s5 += x[i] * x[i] * x[i] * x[i];
s6 += x[i] * y[i];
s7 += x[i] * x[i] * y[i];
}
5. Fitting a parabola
a[0][0] = (float) n;
a[0][1] = s1;
a[0][2] = s3;
a[0][3] = s2;
a[1][0] = s1;
a[1][1] = s3;
a[1][2] = s4;
a[1][3] = s6;
a[2][0] = s3;
a[2][1] = s4;
a[2][2] = s5;
a[2][3] = s7;
guass(a);
}
//Solution using Guass elimination
method
void guass(float a[5][5])
{
float x[5], ratio;
int n = 3, i, j, k;
for(k = 0; k < n - 1; k++)
for(i = k + 1; i < n; i++)
{
ratio = a[i][k] / a[k][k];
for(j = 0; j < n + 1; j++)
a[i][j] -= ratio * a[k][j];
}
5. Fitting a parabola
//Back substitution
x[n-1] = a[n - 1][n] / a[n - 1][n-1];
for(k = n - 2; k >= 0; k--)
{
x[k] = a[k][n];
for(j = k + 1;j < n; j++)
x[k] -= a[k][j] * x[j];
x[k] /= a[k][k];
}
cout <<"The equation for the parabola is: " <<x[2] <<"x^2 + " <<x[1] <<"x +
" <<x[0];
return;
}
6. Trapezoidal rule
void main()
{
float a, b, h, x0, xn, xi;
double sum;
int n, i;
double funct(float a);
cout <<"\n Enter the limits. \n";
cin >>a >>b;
cout <<"Enter the number of
intervals. \n";
cin >>n;
h = (b - a) / n;
x0 = a;
xn = b;
sum = funct (x0) + funct (xn);
cout <<"f(x0) =" <<funct (x0) <<"\n";
for (i = 1; i <= n - 1; i++)
{
xi = a + (i * h);
cout <<"f(x" <<i <<") =" <<funct
(xi) <<"\n";
sum += 2.0 * funct (xi);
}
cout <<"f(x" <<n <<") =" <<funct (xn)
<<"\n";
sum = (h / 2) * sum;
cout <<"The integral is: " <<sum;
}
double funct(float a)
{
double b;
b = exp (-a * a / 2);
return (b);
}
7. Simpsons rule
void main()
{
float a, b, h, x0, xn, xi;
double sum;
int n, i;
double funct(float a);
cout <<"f(x0) = " <<funct (x0) <<"\n";
for (i = 1; i <= n - 1; i++)
{
xi = a + (i * h);
cout <<"\n Enter the limits. \n";
cin >>a >>b;
cout <<"Enter the number of intervals. \n";
cin >>n;
cout <<"f(x" <<i <<") = " <<funct (xi)
<<"\n";
if(i % 2 == 0) sum += 2.0 * funct (xi);
h = (b - a) / n;
x0 = a;
xn = b;
sum = funct (x0) + funct (xn);
else sum += 4.0 * funct (xi);
}
cout <<"f(x" <<n <<") = " <<funct (xn) <<"\n";
sum = (h / 3) * sum;
cout <<"The integral is: " <<sum;
}
7. Simpsons rule
double funct(float a)
{
double b;
b = exp (-a * a / 2);
return (b);
}
8. Gauss elimination method
void main()
{
float a[20][20], ratio, x[20];
int i, j, k, n;
for(j = 0; j < n + 1; j++)
a[i][j] -= ratio * a[k][j];
}
//Back substitution
x[n-1] = a[n - 1][n] / a[n - 1][n-1];
for(k = n - 2; k >= 0; k--)
{
x[k] = a[k][n];
for(j = k + 1;j < n; j++)
x[k] -= a[k][j] * x[j];
x[k] /= a[k][k];
}
for(i = 0; i < n; i++)
cout <<"\nx(" <<i + 1 <<") = " <<x[i];
cout <<"\nEnter the order of the matrix.
\n";
cin >>n;
cout <<"Enter the coefficients & RHS. \n";
//Reading the augmented matrix
for(i = 0; i < n; i++)
for(j = 0; j < n + 1; j++)
cin >>a[i][j];
//Triangularisation
for(k = 0; k < n - 1; k++)
for(i = k + 1; i < n; i++)
{
ratio = a[i][k] / a[k][k];
9. Gauss quadrature
void main()
{
double c[10], x[10];
double sum = 0;
int n, i;
double funct(double a);
cout <<"Enter the number of
Guass points. \n";
cin >>n;
switch(n)
{
case 2:
c[0] = c[1] = 1;
x[0] = -0.57735;
x[1] = 0.57735;
break;
case 3 : c[0] = c[2] = 0.55556;
c[1] = 0.88889;
x[0] = -0.77460;
x[1] = 0;
x[2] = 0.77460;
break;
case 4 :c[0] = c[3] = 0.34785;
c[1] = c[2] = 0.65215;
x[0] = -0.86114;
x[1] = -0.33998;
x[2] = 0.33998;
x[3] = 0.86114;
break;
9. Gauss quadrature
case 5
:c[0] = c[4] = 0.23693;
for(i = 0; i < n; i++)
sum += c[i] * funct(x[i]);
c[1] = c[3] = 0.47863;
cout <<"The integral is: " <<sum;
c[2] = 0.56889;
x[0] = -0.90618;
x[1] = -0.53847;
x[2] = 0;
x[3] = 0.53847;
x[4] = 0.90618;
break;
}
}
double funct(double a)
//Transformed function
{
double b;
b = 3.5 / (3.5 * a + 8.5);
return (b);
}
10. Gauss- Siedel Iteration
void main()
{
float a[10][10],x[10],y[10];
int i,j,n;
clrscr();
/* inputting the coefficients*/
cout<<"enter the number of variables\n";
cin>>n;
cout<<"enter the coefficients\n";
for(i=0;i<n;i++)
{
cout<<"enter row[%d]
elements\n",i+1;
for(j=0;j<=n;j++)
{
cin>>a[i][j];
}
}
/* initialisation of coeficients */
for(i=0;i<n;i++)
x[i]=0;
/* calculation of coefficients */
do
{
for(i=0;i<n;i++)
{ y[i]=x[i];
x[i]=a[i][n];
for(j=0;j<n;j++)
if(i!=j)
x[i]=x[i]-a[i][j]*x[j];
x[i]=x[i]/a[i][i];
}
}while(check(n,x,y));
10. Gauss- Siedel Iteration
/* display of coefficients */
cout<<"The coefficients are:- \n";
for(i=0;i<n;i++)
{
cout<<i<<" "<<x[i];
}
getch();
}
int check(int n, float x[], float y[])
{
int i, flag=0, j;
for(i=0;i<n;i++)
{
if(fabs(x[i]-y[i])>0.00001)
{
flag=1;break;
}
else
flag=0;
}
if (flag==1)
return 1;
else
return 0;
}
