LAB REPORT FOR NUMERICAL METHODS 3

Trapezoidal rule for given function

#include<iostream.h>
#include<conio.h>
#include<math.h>
float F(float x);
void main()
{ int n, i;
float a, b, h, sum, ict;
float F(float x);
cout<<endl<<"TRAPEZOIDAL RULE:"<<endl;
cout<<"input initial value of x"<<endl;
cin>>a;
cout<<"input final value of x:"<<endl;
cin>>b;
cout<<"input the segment width:"<<endl;
cin>>h;
if(a!=b)
{ n= (b-a)/h;
sum = (F(a) + F(b))/2.0;
for(i=1; i<=n-1; i++)
{ sum = sum + F(a+i*h);
}
ict = sum * h;
cout<<endl;
cout<<"Integration between "<<a<<" and "<<b<<" when h= "<<h<<" is "<<ict<<endl;
cout<<endl;
}
else
{ cout<<"exiting"<<endl;
}
getch();
clrscr();
}
float F(float x)
{ float f;
f = 1.0-exp(-x/2.0);
return (f);
}


Trapezoidal rule for tabulated function

#include<iostream.h>
#include<conio.h>
#include<math.h>
#define MAX 15
void main()
{ int n, n1, n2, i;
float a, b, h, sum, ict, x[MAX], y[MAX];
clrscr();
cout<<endl;
cout<<"input the no. of data points:"<<endl;
cin>>n;
cout<<"input the tablated values"<<endl;
cout<<"one set at each line:"<<endl;
for(i=1; i<=n; i++)
cin>>x[i]>>y[i];
cout<<"input initial value of x:"<<endl;
cin>>a;
cout<<"input final value of x:"<<endl;
cin>>b;
cout<<"input width of segments:"<<endl;
cin>>h;
n1=(int)(fabs(a-x[1])/h+1.5);
n2=(int)(fabs(b-x[1])/h+1.5);
sum=0.0;
for(i=n1; i<=n2-1; i++)
sum=sum+y[i]+y[i+1];
ict=sum*h/2.0;
cout<<"integral from "<<a<<" to "<<b<<" is "<<ict<<endl;
getch();
}


Simpson’s 1 by 3 using tabulated value

#include<iostream.h>
#include<conio.h>
#include<math.h>
#include<iomanip.h>
void main()
{ clrscr();
int i, n, n1, n2, m, l;
float x[15], y[15], a, b, h, sum, ics, i1, i2;
cout<<endl;
cout<<"input no. of data points:"<<endl;
cin>>n;
cout<<"input table values set by set"<<endl;
cout<<"one set on each line"<<endl;
for(i=1; i<=n; i++)
cin>>x[i]>>y[i];
cout<<"input initial value of x:"<<endl;
cin>>a;
cout<<"input final value of x:"<<endl;
cin>>b;
cout<<"input the segment width:"<<endl;
cin>>h;
n1 = (int)(fabs(a-x[1])/h+1.5);
n2 = (int)(fabs(b-x[1])/h+1.5);
m=n2-n1;
if(m%2==0)
{ i2=0.0;
l=n2-2;
}
else
{ i2=(y[n2-1]+y[n2])*h/2.0;
l=n2-3;
}
sum=0.0;
for(i=n1; i<=l; i=i+2)
sum=sum+y[i]+4*y[i+1]+y[i+2];
i1=sum*h/3.0;
ics=i1+i2;
cout<<"integral from "<<a<<" to "<<b<<" is "<<ics;
getch();
clrscr();
}

Simpson's 3 by 8 rule integrating a function

#include<iostream.h>
#include<conio.h>
#include<math.h>
#define F(x) 1-exp(-(x)/2.0)
void main()
{ clrscr();
int n, m, i;
float a, b, h, sum, ics, x, f1, f2, f3, f4;
cout<<"Initial value of x"<<endl;
cin>>a;
cout<<"Final value of x"<<endl;
cin>>b;
cout<<"Number of segments n (divisible by 3)"<<endl;
cin>>n;
h=(b-a)/n;
m=n/3;
sum=0.0;
x=a;
f1=F(x);
for(i=0;i<=m-1;i++)
{ f2=F((x+h)); c
f3=F((x+2*h));
f4=F((x+3*h));
sum=sum+f1+3*f2+3*f3+f4;
f1=f4;
x=x+3*h;
}
ics=sum*3*h/8.0;
cout<<"Integral from "<<a<<" to "<<b<<endl;
cout<<"when h= "<<h<<" is "<<ics<<endl;
getch();
}

LAB REPORT FOR NUMERICAL METHODS 4

Gauss siedel iteration method

#include<iostream.h>
#include<conio.h>
#include<iomanip.h>
#include<math.h>
#define MAXIT 50
#define EPS 0.000001
void gaseid(int n, float a[10][10], float b[10], float x[10], int *count, int *status);
void main()
{ float a[10][10], b[10], x[10];
int i, j, n, count, status;
cout<<"** SOLUTION BY GUASS SEIDEL ITERATION METHOD **"<<endl;
cout<<"input the size of the system:"<<endl;
cin>>n;
cout<<"input coefficients, a(i.j)"<<endl;
cout<<"one row on each line"<<endl;
for(i=1; i<=n; i++)
for(j=1; j<=n; j++)
cin>>a[i][j];
cout<<"input vector b:"<<endl;
for(i=1; i<=n; i++)
cin>>b[i];
gaseid(n, a, b, x, &count, &status);
if(status==2)
{ cout<<"no CONVERGENCE in "<<MAXIT<<" iterations."<<endl<<endl<<endl;
}
else
{ cout<<"SOLUTION VECTOR X"<<endl;
for(i=1; i<=n; i++)
cout<<setw(15.6)<<x[i]<<endl;
cout<<"iterations= "<<count;
}
getch();
}
void gaseid(int n, float a[10][10], float b[10], float x[10], int *count, int *status)
{ int i, j, key;
float sum, x0[10];
for(i=1; i<=n; i++)
x0[i]=b[i]/a[i][i];
*count=1;
begin:
key=0;
for(i=1; i<=n; i++)
{ sum=b[i];
for(j=1; j<=n; j++)
{ if(i==j)
continue;
sum=sum-a[i][j]*x0[j];
}
x[i]=sum/a[i][i];
if(key==0)
{ if(fabs((x[i]-x0[i])/x[i])>EPS)
key=1;
}
}
if(key==1)
{ if(*count==MAXIT)
{ *status=2;
return;
}
else
{ *status=1;
for(i=1; i<=n; i++)
x0[i]=x[i];
}
*count=*count+1;
goto begin;
}
return;
}



Euler’s Method
#include<iostream.h>
#include<conio.h>
#include<math.h>
#include<iomanip.h>
float func(float x, float y);
void main()
{ clrscr();
int i, n;
float x, y, xp, h, dy;
float func(float, float);
cout<<endl;
cout<<"SOLUTION BY THE EULER's METHOD"<<endl;
cout<<"input the initial values of x and y:"<<endl;
cin>>x>>y;
cout<<"input the x at which y is required:"<<endl;
cin>>xp;
cout<<"input the step size, h:"<<endl;
cin>>h;
n = (int) ((xp-x)/h + 0.5);
for(i=1; i<=n; i++)
{ dy = h * func(x, y);
x = x + h;
y = y + dy;
cout<<setw(5)<<i<<setw(10.6)<<x<<setw(10.6)<<y<<endl;
}
cout<<"the value of y at x = "<<x<<" is "<<y<<endl;
getch();
clrscr();
}
float func(float x, float y)
{ float f;
f = x + y + x*y;
return (f);
}




Heun’s Method
#include<iostream.h>
#include<conio.h>
#include<math.h>
#include<iomanip.h>
float func(float x, float y);
void main()
{ clrscr();
int i, n;
float x, y, xp, h, m1, m2;
float func(float, float);
cout<<endl;
cout<<"SOLUTION BY THE HEUN's METHOD"<<endl;
cout<<endl;
cout<<"input initial values of x and y:"<<endl;
cin>>x>>y;
cout<<"input x at which y is required:"<<endl;
cin>>xp;
cout<<"input the step size, h:"<<endl;
cin>>h;
n = (int) ((xp - x)/h + 0.5);
for(i=1; i<=n; i++)
{ m1 = func(x, y);
m2 = func(x+h, y+m1*h);
x = x + h;
y = y + 0.5 * h * (m1 + m2);
cout<<setw(7)<<i<<setw(15.9)<<x<<setw(15.9)<<y<<endl;
}
cout<<"the value of y at x = "<<x<<" is "<<y;
getch();
clrscr();
}
float func(float x, float y)
{ float f;
f = -y/(2*y + 1);
return(f);
}


Runge-Kutta Method
#include<iostream.h>
#include<conio.h>
#include<math.h>
#include<iomanip.h>
float func(float x, float y);
void main()
{ clrscr();
int i, n;
float x, y, xp, h, m1, m2, m3, m4;
float func(float, float);
cout<<endl;
cout<<"******************************************"<<endl;
cout<<"SOLUTION BY 4th ORDER RUNGE- KUTTA METHOD:"<<endl;
cout<<"******************************************"<<endl;
cout<<"input initial values of x and y:"<<endl;
cin>>x>>y;
cout<<"input x at which y is required:"<<endl;
cin>>xp;
cout<<"input step size, h:"<<endl;
cin>>h;
n = (int) ((xp - x)/h + 0.5);
cout<<endl;
cout<<"------------------------------------"<<endl;
cout<<setw(5)<<"STEP"<<setw(15.9)<<"X"<<setw(15.9)<<"Y"<<endl;
cout<<"------------------------------------"<<endl;
for(i=1; i<=n; i++)
{ m1 = func(x, y);
m2 = func(x + 0.5*h, y + 0.5*m1*h);
m3 = func(x + 0.5*h, y + 0.5*m2*h);
m4 = func(x+h, y + m3*h);
x = x + h;
y = y + (m1 + 2.0*m2 + 2.0*m3 + m4) * h/6.0;
cout<<setw(5)<<i<<setw(15.9)<<x<<setw(15.9)<<y<<endl;
}
cout<<"the value of y at x = "<<x<<" is "<<y<<endl;
getch();
clrscr();
}
float func(float x, float y)
{ float f;
f = y + sqrt(y);
return(f);
}

LAB REPORT FOR NUMERICAL METHODS 2

Lagrange's Interpolation

#include<iostream.h>
#include<math.h>
#include<conio.h>
#define MAX 10
void main()
{ clrscr();
int n, i, j;
float x[MAX], f[MAX], fp, lf, sum, xp;
cout<<"input no. of data points:"<<endl;
cin>>n;
cout<<"input values of x and fx"<<endl;
cout<<"one set on each line"<<endl;
for(i=1; i<=n; i++)
cin>>x[i]>>f[i];
cout<<"input x where the interpolation is required:"<<endl;
cin>>xp;
sum=0.0;
for(i=1; i<=n; i++)
{ lf=1.0;
for(j=1; j<=n; j++)
{ if(i!=j)
lf=lf*(xp-x[j])/(x[i]-x[j]);
}
sum= sum+lf*f[i];
}
fp= sum;
cout<<" LAGRANGIAN INTERPOLATION "<<endl;
cout<<"interpolated function value:"<<endl;
cout<<"at x= "<<xp<<" it is "<<fp<<endl;
getch();
}


Newton’s Interpolation

#include<iostream.h>
#include<conio.h>
void main()
{ clrscr();
int i, j, n;
float xp, fp, sum, pi, x[10], f[10], a[10], d[10][10];
cout<<"input numbers of data points:"<<endl;
cin>>n;
cout<<"Input values of x and fx"<<endl;
cout<<"one set on each line"<<endl;
for(i=1; i<=n; i++)
cin>>x[i]>>f[i];
for(i=1; i<=n; i++)
d[i][1]=f[i];
for(j=2; j<=n; j++)
for(i=1; i<=n-j+1; i++)
d[i][j]=(d[i+1][j-1]-d[i][j-1])/(x[i+j-1]-x[i]);
for(j=1; j<=n; j++)
a[j]=d[1][j];
cout<<"enter xp where interpolation is required:"<<endl;
cin>>xp; csitnepal
Source: www.csitnepal.com Page 19
sum=a[1];
for(i=2; i<=n; i++)
{ pi =1.0;
for(j=1; j<=i-1; j++)
pi=pi*(xp-x[j]);
sum=sum+a[i]*pi;
}
fp=sum;
cout<<"NEWTON's INTERPOLATION"<<endl;
cout<<"***********************"<<endl;
cout<<endl;
cout<<"interpolation values"<<endl;
cout<<"at x="<<xp<<" it is"<<fp<<endl;
getch();
}

Fitting a straight line

#include<iostream.h>
#include<conio.h>
#include<math.h>
#define MAX 10
#define EPS 0.000001
void main()
{ clrscr();
int i, n;
float x[10], y[10];
float sumx, sumy, sumxx, sumxy, xmean, ymean;
float denom;
float a,b;
cout<<"LINEAR REGRESSION METHOD"<<endl;
cout<<endl;
cout<<"input the no. of data points:"<<endl;
cin>>n;
cout<<"input x and y values"<<endl;
cout<<"one set on each line:"<<endl;
for(i=1; i<=n; i++)
cin>>x[i]>>y[i];
sumx=0.0; csitnepal
Source: www.csitnepal.com Page 22
sumy=0.0;
sumxx=0.0;
sumxy=0.0;
for(i=1; i<=n; i++)
{ sumx=sumx+ x[i];
sumy=sumy+ y[i];
sumxx=sumxx+ x[i]*x[i];
sumxy= sumxy+ x[i]*y[i];
}
xmean=sumx/n;
ymean=sumy/n;
denom=n*sumxx-sumx*sumx;
if(fabs(denom)>EPS)
{ b=(n*sumxy-sumx*sumy)/denom;
a=ymean-b*xmean;
cout<<"linear regression line y=a+bx"<<endl;
cout<<"the coefficients are:"<<endl;
cout<<"a= "<<a<<endl;
cout<<"b= "<<b<<endl;
}
else
cout<<"NO SOLUTION"<<endl;
getch();
}

Differentiation using newton's interpolating polynomial

#include<iostream.h>
#include<conio.h>
#include<math.h>
void main()
{ int i, j, k, n;
float x[10], f[10], a[10], d[10][10], xp, dif, sum, p;
clrscr();
cout<<"NUMERICAL DIFFERENTIATION USING NEWTON's POLYNOMIAL"<<endl;
cout<<"input no. of data points:"<<endl;
cin>>n;
cout<<"input values of x and f(x)"<<endl;
cout<<"one set on each line"<<endl;
for(i=1; i<=n; i++)
cin>>x[i]>>f[i];
for(i=1; i<=n; i++)
d[i][1] = f[i];
for(j=2; j<=n; j++)
for(i=1; i<=n-j+1; i++)
d[i][j]=(d[i+1][j-1]-d[i][j-1])/(x[i+j-1]-x[i]);
for(j=1; j<=n; j++)
a[j]=d[1][j];
cout<<"input xp where derivative is required:"<<endl;
cin>>xp;
dif =a[2];
for(k=3; k<=n; k++)
{ sum= 0.0;
for(i=1; i<=k-1; i++)
{ p= 1.0;
for(j=1; j<=k-1; j++)
{ if (i==j)
continue;
p= p*(xp-x[j]);
}
sum= sum+p;
}
dif=dif+a[k]*sum;
}
cout<<"derivative at x= "<<xp<<" is "<<dif<<endl;
getch();
}

LAB REPORT FOR NUMERICAL METHODS 1

Bisection Method

#include<iostream.h>
#include<conio.h>
#include<math.h>
#define EPS 0.000001
#define F(x) log(x)-cos(x)
void bim(float*a, float*b, float*root, float*s, int*count);
void main()
{ int count;
float a, b, root, s;
cout<<"SOLUTION BY BISECTION METHOD:"<<endl;
cout<<"input starting values:"<<endl;
cin>>a>>b;
bim(&a, &b, &root, &s, &count);
if (s==0)
{ cout<<"starting points donot bracket any root"<<endl;
cout<<"check whether they bracket EVEN roots."<<endl;
}
else
{ cout<<"root="<<root<<endl;
float fun = F(root);
cout<<"f(root)="<<fun<<endl;
cout<<"Iterations="<<count<<endl;
}
getch();
clrscr();
}
void bim(float*a, float*b, float*root, float*s, int*count)
{ float x1, x2, x0, f0, f1, f2;
x1=*a;
x2=*b;
f1=F(x1);
f2=F(x2);
if(f1*f2> 0)
{ *s=0;
return;
}
else
{ *count=0;
begin:
x0=(x1+x2)/2.0;
f0=F(x0);
if(f0==0)
{ *s=1;
*root=x0;
return;
}
if(f1*f0<0)
{ x2=x0;
}
else
{ x1=x0;
f1=f0;
}
if(fabs((x2-x1)/x2) < EPS)
{ *s=1;
*root=(x1+x2)/2.0;
return;
}
else
{ *count=*count+1;
goto begin;
}
}
}




Secant Method

#include<iostream.h>
#include<conio.h>
#include<math.h>
#define EPS 0.000001
#define MAXIT 50
#define F(x) exp(x)-x-2
int sec(float*a, float*b, float*x1, float*x2, float*root, int*count, int*status);
void main()
{ float a, b, root, x1, x2, fr;
int count, status;
clrscr();
cout<<"SOLUTION BY SECANT METHOD"<<endl;
cout<<"Input two starting points:"<<endl;
cin>>a>>b;
sec(&a, &b, &x1, &x2, &root, &count, &status);
if(status==1)
{ cout<<"DIVISION BY ZERO"<<endl;
cout<<"last x1="<<x1<<endl<<"last x2="<<x2<<endl;
cout<<"number of iterations="<<count<<endl;
}
if(status==2) csitnepal
Source: www.csitnepal.com Page 6
{ cout<<"NO CONVERGENCE IN"<<MAXIT<<" ITERATIONS."<<endl;
}
else
{ cout<<"root="<<root<<endl;
fr=F(root);
cout<<"function Value at root="<<fr<<endl;
cout<<"no. of iterations="<<count<<endl;
}
getch();
}
int sec(float*a, float*b, float*x1, float*x2, float*root, int*count, int*status)
{ float x3, f1, f2, error;
*x1=*a;
*x2=*b;
f1=F(*a);
f2=F(*b);
*count=1;
begin:
if (fabs(f1-f2)<=1.E-10)
{ *status=1;
return 0;
}
x3=*x2-f2*(*x2-*x1)/(f2-f1);
error=fabs((x3-*x2)/x3);
if (error>EPS)
{ if (*count==MAXIT)
{ *status=2;
return 0;
}
else
{ *x1=*x2;
}
*x2=x3;
f1=f2;
f2=F(x3);
*count=*count+1;
goto begin;
}
else
{ *root=x3;
*status=3;
return 0;
}
}



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<<"SOLUTION BY NEWTON RAPHSON'S METHOD"<<endl;
cout<<"input initial value of x:"<<endl;
cin>>x0;
count=1;
begin:
fx=F(x0);
fdx=FD(x0);
xn=x0-fx/fdx;
if (fabs((xn-x0)/xn) < EPS)
{ cout<<"root="<<xn<<endl;
fxn =F(xn);
cout<<"function value="<<fxn<<endl;
cout<<"no. of iterations="<<count<<endl;
}
else
{ x0=xn;
count=count+1;
if(count<MAXIT)
{ goto begin;
}
else
{ cout<<"SOLUTION DOESNOT CONVERGE."<<endl;
cout<<"iterations="<<MAXIT<<endl;
}
}
getch();
}
csitnepal



Fixed-point Iteration Method

#include<iostream.h>
#include<conio.h>
#include<iomanip.h>
#include<math.h>
#define EPS 0.000001
#define G(x) 2.0-(x)*(x)
void main()
{ clrscr();
int MAXIT, i;
float x0, x, error;
cout<<"SOLUTON BY FIXED POINT METHOD"<<endl;
cout<<"input initial estimate of a root:"<<endl;
cin>>x0;
cout<<"Maximum iterations allowed:"<<endl;
cin>>MAXIT;
cout<<"iteration value of X error:"<<endl;
for(i=1; i<=MAXIT; i++)
{ x = G(x0);
error=fabs((x-x0)/x);
cout<<setw(10)<<i<<setw(10)<<x<<setw(10.8)<<error<<endl;
if(error < EPS)
goto end;
else
x0=x;
}
cout<<"process doesnot converge to a root"<<endl;
cout<<"exit from the iteratoion loop"<<endl;
end:
;
getch();
}