Problem: 1 1.

(a) Program to solve nonlinear algebraic or transcendental function with example

by Bisection method.

Example: Find a real root of the equation x3 – 2x − 5 = 0.

Solution:
Let ƒ (x) = x3 – 2x − 5.
Then
ƒ (2) = −1 and ƒ (3) = 16.
Hence a root lies between 2 and 3 and we take
x1=(2+3)/2=2.5
Since ƒ (X1) = ƒ (2.5) = 5.6250, the root lies between 2 and 2.25. Hence
X2=(2+2.5)/2=2.25
Now, ƒ (X2) = 1.890625, the root lies between 2 and 2.25.
Therefore,
X3=(2+2.25)/2=2.125
Since ƒ (X3) = 0.3457, the root lies between 2 and 2.125.
Therefore,
X4=(2+2.125)/2=2.0625
Proceeding in this way, we obtain the successive approzimations:
X5 = 2.09375, x6 = 2.10938, x7 = 2.10156,
X8= 2.09766, x9 = 2.09570, x10 = 2.09473,
X11 = 2.09424, …

Hence the root correct upto three decimal places is 2.09424

Program:

#include<bits/stdc++.h>
using namespace std;
double f(double x)
{
double aa=pow(x,3)-2.0*x-5.0;
return aa;
}
int main()
{
start:
cout.precision(7);
cout.setf(ios::fixed);
double a,b,c,ep;
int i;
k:cout<<"Enter the value of interval\n";
cin>>a;
cin>>b;
cout<<"Enter the value of desired acuracy\n";
cin>>ep;
if(f(a)*f(b)>0)

1
 {
cout<<"Please enter another interval\n";
goto k;
}
else
{
while(abs(a-b)>=ep)
{
c=(a+b)/2.0;
if(f(c)==0)
{
cout<<"The root is="<<c<<endl;
break;
}
if(f(c)*f(a)>0)
{
a=c;
}
else
{
b=c;
}
i++;
}
cout<<"\nThe root is:"<<c<<endl;
cout<<"\nThe number of iteration is:"<<i<<endl;
}
cout << "Press (0) to exit or (1) to continue the next test: ";
int ch; cin >> ch;
if(ch==1) {goto start;}
}

Output:
Enter the value of interval
2 3
Enter the value of desired acuracy
0.001

The root is : 2.0947266

The number of iteration is:10

Comments: Since the root obtained by the code is almost the same as the root that we found
earliest by manual calculation. So this code could be implemented to solve other algebraic or
transcendental functions to find the root using the bisection method.

Problem: 2 1.(b) Program to solve nonlinear algebraic or transcendental function with example by
iterative method.

Example: Find a real root, correct to three decimal places, of the equation 2x − 3 = cos x.

Solution:

2
 We rewrite the given equation as
x=0.5(cos x+3)
Choosing x0 = 1.5, we obtain successively:
x1 = 0.5 (cos 1.5 +3) = 1.5354
x2 = 0.5 (cos 1.5354 +3)= 1.5177
x3 = 0.5 (cos 1.5177 +3) = 1.5265
x4 = 0.5 (cos 1.5265 + 3) = 1.5221
x5 = 0.5 (cos 1.5221 +3) = 1.5243
x6 = 0.5 (cos 1.5243 +3) = 1.5232
x7 = 0.5 (cos 1.5232 +3) = 1.5238
Hence the root correct upto three decimal places is 1.524.

Program:

#include<bits/stdc++.h>
using namespace std;
#define f(x) (0.5*(cos(x)+3))
int main()
{
double x0, x, ac,r;
int i=1;

cout<< setprecision(6)<< fixed;

cout<<"Enter initial guess from (3/2,pi/2) interval \n";
cin>>x0;
cout<<"Please enter the desired accuracy\n ";
cin>>ac;
do
{
x = f(x0);
r=x0;
x0= x;
cout<< endl<<"x["<<i<<"] ="<< x;
i++;
}while( fabs(x-r)> ac);
cout<<"The root is :"<<x<<endl;
return(0);
}

Output:
Enter initial guess from (3/2,pi/2) interval
1.5
Please enter the desired accuracy
0.001

x[1] =1.535369
x[2] =1.517710
x[3] =1.526531

3
 x[4] =1.522126
x[5] =1.524326
x[6] =1.523227
x[7] =1.523776
The root is :1.523776.

Comments: Since the root obtained by the code is almost the same as the root that we found
earliest by manual calculation. So this code could be implemented to solve other algebraic or
transcendental functions to find the root using the iterative method.

Problem: 3 1.(c) Program to solve nonlinear algebraic or transcendental function with example by
Newton-Raphson method.

Example: Using Newton−Raphson method, find a real root, correct to 3 decimal

places, of the equation sin x= x/2 given that the root lies between π/2 and π.

Solution: Let f(x) = sin x − x/2

Then f ‘(x) = cos x − 1/2

Choosing x0 = π/2,we obtain

X1 = π/2 − ( sin π/2− π/4)/ −0.5 =2,
X2 =2 − ( sin 2− 1)/(cos 2 − 0.5) = 1.9010,
X3 =1.9010 − ( sin 1.9010 – 1)/( cos 1.9010 – 0.5) = 1.8955 ,
Similarly x4 = 1.8954 , x5 = 1.8955 , …….
Hence the required root is x = 1.896 .

Program:

#include<bits/stdc++.h>
using namespace std;

double f(double x)
{
double phi=sin(x)-x/2.0;
return phi;
}
double g(double x)
{
double deri=cos(x)-1.0/2.0;
return deri;
}
int main()
{
cout.precision(7);
cout.setf(ios::fixed);
double x,accuracy;
int i;
cout<<"please enter the initial approximation\n";
cin>>x;
cout<<"please enter the required accuray\n";
cin>>accuracy;

4
 if(f(x)==0)
{
cout<<"the root is ="<<x;
}
else
{
i=1;
do
{
x=x-f(x)/g(x);

cout<<x<<endl;
i++;
}while(abs(f(x)-0)>=accuracy);
}
cout<<"the root is:"<<x<<endl;
cout<<"the number of iteration is:"<<i<<endl;
}
Output:
please enter the initial approximation
1.5708
please enter the required accuray
0.001

1.9999968
1.9009953
1.8955116
the root is:1.8955116
the number of iteration is:4

Comments: It's great to hear that your code is producing results that are close to the manually
calculated roots. The Newton-Raphson method is indeed a powerful numerical technique for finding
roots of algebraic or transcendental functions. As long as you have a function that is continuous and
differentiable in the neighborhood of the root you want to find, you can implement the Newton-
Raphson method to approximate that root.

Problem: 4 1.(d) Program to solve a system of nonlinear equations (two variables only) with
example by Newton-Raphson method.

Example : Solve the system

sin x − y = −0.9793

cos y − x = − 0.6703
with x0 = 0.5 and y0 = 1.5 as the initial approzimation.
Solution :
We have
ƒ (x, y) = sin x − y + 0.9793
and g(x, y) = cos y − x + 0.6703
For first iteration, we have
f0 = - 0.0413 , g 0 = 0.3410
D = fx g y − fy g x =cos (0.5) ( − sin 1.5 ) – 1 = − 1.8754

5
 fgy −gfx gfx −𝑓gx
h= D
= − 0.1505 k= D
= − 0.0908

therefore
x = 0.5 + 0.1505 = 0.6505
y = 1.5 + 0.0908 = 1.5908
For the second iteration, we have x0 = 0.6505 and y0 = 1.5908. Then we
obtain
D = −1.7956
Also h = −0.003181 and k = 0.003384.
Hence the new approzimation is
x = 0.6505 + 0.0032 = 0.6537,
y = 1.5908 − 0.0034 = 1.5874.

Program:

#include <iostream>
#include <cmath>
#include <iomanip>
using namespace std;
#define f(x,y) (sin(x)-y+0.9793)
#define g(x,y) (cos(y)-x+0.6703)
#define dfx(x,y) (cos(x))
#define dfy(x,y) (-1)
#define dgx(x,y) (-1)
#define dgy(x,y) (-sin(y))
int main()
{
cout.precision(6);
cout.setf(ios::fixed);
double x, x0, y, y0, d, h, k, ac;
cout<<"Enter the initial guesses x0 and y0 respectively:\n";
m:
cin>>x0>>y0;
cout<<"Enter the desired accuracy:\n"<<endl;
cin>>ac;
do
{
x=x0,y=y0;
d=dfx(x,y)*dgy(x,y)-dgx(x,y)*dfy(x,y);
if (d==0)
{
cout<<"The system doesn't converge for this initial guess. Please, give new
guess:"<<endl;
goto m;
}
else
{

6
 h=(g(x,y)*dfy(x,y)-f(x,y)*dgy(x,y))/d;
k=(f(x,y)*dgx(x,y)-g(x,y)*dfx(x,y))/d;
x0=x+h, y0=y+k;
}
}while(fabs(x-x0)>ac || fabs(y-y0)>ac);
cout<<"Root of the equation: x= "<<x0<<" y= "<<y0<<endl;
return 0;
}

Output:
Enter the initial guesses x0 and y0 respectively:
0.5 1.5
Enter the desired accuracy:
0.001
Root of the equation: x= 0.653683 y= 1.587414

Comments: The Newton-Raphson method is a powerful numerical technique for finding roots of
algebraic or transcendental functions. It can also be extended to two-variable functions. To
implement the Newton-Raphson method for two-variable functions. It's great to hear that your code
is producing results that are close to the manually calculated roots. The Newton-Raphson method is
indeed a powerful numerical technique for finding roots of algebraic or transcendental functions. As
long as you have a function that is continuous and differentiable in the neighborhood of the root you
want to find, you can implement the Newton-Raphson method to approximate that root.

Problem: 5 2.(a) Program to calculate Interpolation by Newton’s backward formula with

example.

Example : Values of x (in degrees) and sin x are given in the following table:

x (in degrees) Sin x

15 0.2588190
20 0.3420201
25 0.4226183
30 0.5
35 0.5735764
40 0.6427876
Determine the value of sin 380 .
Solution :
The difference table is :

x Sin x ∆ ∆2 ∆3 ∆4 ∆5
15 0.2588190
20 0.3420201 0.0832011
25 0.4226183 0.0805982 -0.0026029
30 0.5 0.0773817 -0.0032165 -0.0006136
35 0.5735764 0.0735764 -0.0038053 -0.0005888 0.0000248
40 0.6427876 0.0692112 -0.0043652 -0.0005599 0.0000289 0.0000041

To find sin 38°, we use Newton’s backward difference formula with xn = 40 and x = 38. This
gives,
x−x 38−40
p = h n = 5 = −0.4

7
 Hence we obtain
(−0.4)(−0.4+1)
𝑦(38) = 0.6427876 − 0.4(0.069212) + (−0.0043652)
2

(−0.4)(−0.4 + 1)(0.4 + 2)
+ (−0.000559)
6
(−0.4)(−0.4 + 1)(0.4 + 2)(−0.4 + 3)
+ (0.0000289)
24
(−0.4)(−0.4 + 1)(0.4 + 2)(−0.4 + 3)(−0.4 + 4)
+ (0.0000041)
120
=0.6427876 − 0.02768448 + 0.00052382 + 0.00003583 − 0.00000120

=0.6156614

Program :

#include<iostream>
#include<iomanip>
using namespace std;
int main()
{
cout<<setprecision(5)<<fixed;

int i=0,j=0,n,k=0;
cout<<"\nEnter the number of values to be entered.\n";
cin>>n;
double x[n], y[n][n];
cout<<"\nEnter the values of x\n";
for (i=0;i<n;i++)
cin>>x[i];
cout<<"\nEnter the values of y\n";
for (i=0;i<n;i++)
cin>>y[i][0];
// Table
for (j=1;j<n;j++)
{
k++;
for (i=k;i<n;i++)
{
y[i][j]=y[i][j-1]-y[i-1][j-1];
}
}
cout<<"\n The Backward Difference Table is as follows: \n\n";
cout<<"x"<<setw(15)<<"y"<<setw(15);
for (i=1;i<n;i++)
cout<<"y"<<i<<setw(15);
cout<<"\n-----------------------------------------------------------------------\n";
for (i=0;i<n;i++)
{
cout<<x[i]<<setw(15);
for (j=0;j<=i;j++)
{

8
 cout<<y[i][j];
cout<<setw(15);
}
cout<<"\n";
}
double h, fx, p, sum, tt=1.0;
cout<<"\nEnter the value of x at which the value of y to be calculated\n";
cin>>fx;
h=x[1]-x[0];
p=(fx-x[n-1])/h;
sum=y[n-1][0];
for(j=1;j<n;j++)
{
tt=tt*(p+j-1)/j;
sum+=tt*y[n-1][j];
}
cout<<"\nThe value of y at x=38 is\n\n";
cout<<"y("<<fx<<")="<<sum;
return 0;
}

Output:

Enter the number of values to be entered.

6
Enter the values of x
15 20 25 30 35 40
Enter the values of y
0.2588190
0.3420201
0.4226183
0.5
0.5735764
0.6427876

The Backward Difference Table is as follows:

x y y1 y2 y3 y4 y5
-----------------------------------------------------------------------
15.00000 0.25882
20.00000 0.34202 0.08320
25.00000 0.42262 0.08060 -0.00260
30.00000 0.50000 0.07738 -0.00322 -0.00061
35.00000 0.57358 0.07358 -0.00381 -0.00059 0.00002
40.00000 0.64279 0.06921 -0.00437 -0.00056 0.00003 0.00000

Enter the value of x at which the value of y to be calculated

38

The value of y at x=38 is

9
Y (38.00000)=0.61566

Comments: Since the obtained value by the code is almost the same as we found earlier by
manual calculation. So this code could be implemented to solve other algebraic or
transcendental functions to find the required value using Newton's backward formula.

Problem: 6 2. (b) Program to calculate Interpolation by Newton’s Forward formula with

example .

Example: The table gives the values of tan x for 0.10 ≤ x ≤ 0. 30 . Find the value of tan 0. 12 .

x Y = tan x
0.10 0.1003
0.15 0.1511
0.20 0.2027
0.25 0.2553
0.30 0.3093

Solution: The difference table is

x y ∆ ∆2 ∆3 ∆4
0.10 0.1003 0.0508 0.0008 0.0002 0.0002
0.15 0.1511 0.0516 0.0010 0.0004
0.20 0.2027 0.0526 0.0014
0.25 0.2553 0.0540
0.30 0.3093

By Newton’s forward difference formula :

0.4 (0.4−1) 0.4 (0.4−1)(0.4−2)

tan(0.12) = 0.1003 + 0.4 (0.0508) + (0.0008) + (0.0002)
2 6
0.4 (0.4−1)(0.4−2)(0.4−3)
+ (0.0002)
24
= 0.1205

10
Program :

#include<iostream>
#include<iomanip>
using namespace std;
int main()
{
cout<<setprecision(5)<<fixed;

int i=0,j=0,n,k=0;
cout<<"\nEnter the number of values to be entered.\n";
cin>>n;
double x[n], y[n][n];
cout<<"\nEnter the values of x\n";
for (i=0;i<n;i++)
cin>>x[i];
cout<<"\nEnter the values of y\n";
for (i=0;i<n;i++)
cin>>y[i][0];

// Table calculation
for (j=1;j<n;j++)
for (i=0;i<n-j;i++)
{
y[i][j]=y[i+1][j-1]-y[i][j-1];
}
cout<<"\n The Forward Difference Table is as follows: \n\n";
cout<<"x"<<setw(15)<<"y"<<setw(15);
for (i=1;i<n;i++)
cout<<"y"<<i<<setw(15);
cout<<"\n-----------------------------------------------------------------------\n";
k=n;
for (i=0;i<n;i++)
{
cout<<x[i]<<setw(15);
for (j=0;j<k;j++)
{
cout<<y[i][j];
cout<<setw(15);
}
cout<<"\n";
k--;
}
double h, fx,p, sum, tt=1.0;
cout<<"\nEnter the value of x at which the value of y to be calculated\n";
cin>>fx;
h=x[1]-x[0];
p=(fx-x[0])/h;
sum=y[0][0];
for(j=1;j<n;j++)
{

11
 tt=tt*(p-j+1)/j;
sum+=tt*y[0][j];
}
cout<<"y("<<fx<<")="<<sum;
return 0;
}

Output:
Enter the values of x
0.10 0.15 0.20 0.25 0.30

Enter the values of y

0.1003
0.1511
0.2027
0.2553
0.3093

The Forward Difference Table is as follows:

x y y1 y2 y3 y4
-----------------------------------------------------------------------
0.10000 0.10030 0.05080 0.00080 0.00020 0.00020
0.15000 0.15110 0.05160 0.00100 0.00040
0.20000 0.20270 0.05260 0.00140
0.25000 0.25530 0.05400
0.30000 0.30930

Enter the value of x at which the value of y to be calculated

0.12
Y (0.12000) = 0.12053

Comments: Since the obtained the required value by the code is almost the same as we found
earlier by manual calculation. So this code could be implemented to solve other algebraic or
transcendental functions to find the required value using Newton's Forward formula.

Problem: 7 2.(c) Program to calculate Interpolation and extrapolation by Lagrange’s formula with
example.

Example: consider thevalues of x and y

x 4 12 19
y 1 3 4

Find y when x= 7.

Solution :
Applying Lagrange interpolation formula

(−5)(−12) (3)(−12) (3)(−5)

y = (−8)(−15) (1) + (8)(−7)
(−3) + (15)(7) (4)

12
 1 27 4
= + +
2 14 7

=1.86

Program :

#include<iostream>
#include<iomanip>
#include<cmath>
using namespace std;
int main()
{
cout.precision(7);
cout.setf(ios::fixed);
int i, j, n;
cout<<"Enter the number of values to be entered\n";
cin>>n;
double x[n], y[n];
cout<<"Enter the of values of x\n";
for(i=0; i<n; i++)
cin>>x[i];
cout<<"Enter the of values of y\n";
for(i=0; i<n; i++)
cin>>y[i];

//code of interpolation starts here

double fx, sum=0, temp;
cout<<"Enter the values of x at which y to be calculated\n";
cin>>fx;
for(j=0; j<n; j++)
{
temp=1;
for(i=0; i<n; i++)
{
if(i!=j)
temp=temp*(fx-x[i])/(x[j]-x[i]);
}
sum=sum+temp*y[j];
}
cout<<"The value of y at x="<<fx<<" is: "<<sum<<endl;
return 0;
}

Output:
Enter the number of values to be entered
3
Enter the of values of x
4 12 19
Enter the of values of y
1 3 4
Enter the values of x at which y to be calculated

13
 7
The value of y at x = 7.0000000 is : 1.8571429

Comments: Since the obtained the required value by the code is almost the same as we found
earlier by manual calculation. So this code could be implemented to solve other algebraic or
transcendental functions to find the required value using Interpolation and extrapolation by
Lagrange’s formula.

Problem: 8 2.(d) Program to calculate Interpolation and extrapolation by Newton’s divided

difference formula with example.

Example: The divided difference table is

x 𝑙𝑜𝑔10 𝑥 [𝑥0 , 𝑥1 ] [𝑥0 , 𝑥1 , 𝑥2 ]

300 2.4771 0.00145 0.00001
304 2.4829 0.00140 0
305 2.4843 0.00140
307 2.4871

Find the value of log10 301 .

Solution:
Newton’s divided difference formula is ,
y = y0 + (x − x0 )[x0 , x1 ] + (x − x0 )(x − x1 )[x0 , x1 , x2 ]

log10 301 = 2.4771 + 0.00145 + (−3)(−0.00001)

= 2.4786

Program:

#include<iostream>
#include<iomanip>
#include<cmath>
using namespace std;
int main()
{
cout.precision(6);
cout.setf(ios::fixed);
int i, j, k, n;
cout<<"Enter the number of values to be entered\n";
cin>>n;
double x[n], y[n][n];
cout<<"Enter the of values of x\n";
for(i=0; i<n; i++)
cin>>x[i];
cout<<"Enter the of values of y\n";
for(i=0; i<n; i++)
cin>>y[i][0];
//Calculate difference table
for(j=1; j<n; j++)
{

14
 for(i=0; i<n-j; i++)
{
y[i][j]=(y[i+1][j-1]-y[i][j-1])/(x[i+j]-x[i]);
}
}
//Print difference table
cout<<"\nThe Divided difference table is as follows:\n";
cout<<"x"<<setw(20)<<"y"<<setw(20);
for(i=1; i<n; i++)
cout<<"d"<<i<<"y"<<setw(20);
cout<<"\n-------------------------------------------------------\n";
k=n;
for(i=0; i<n; i++)
{
cout<<x[i]<<setw(20);
for(j=0; j<k; j++)
{
cout<<y[i][j]<<setw(20);
}
cout<<"\n";
k--;
}

//Code of interpolation
double fx, sum=y[0][0], temp=1.0;
cout<<"Enter the values of x at which y to be calculated\n";
cin>>fx;
for(j=1; j<n; j++)
{
temp=temp*(fx-x[j-1]);
sum=sum+temp*y[0][j];
}
cout<<"The value of y at x="<<fx<<" is: "<<sum<<endl;

double per_err;
per_err=(abs(log10(fx)-sum)/log10(fx))*100;
cout<<"Percent Error = "<<per_err<<" % "<<endl;
return 0;
}

Output:
Enter the number of values to be entered
4
Enter the of values of x
300 304 305 307
Enter the of values of y
2.4771
2.4829
2.4843
2.4871

15
 The Divided difference table is as follows:
x y d1y d2y d3y
-------------------------------------------------------
300.000000 2.477100 0.001450 -0.000010 0.000001
304.000000 2.482900 0.001400 -0.000000
305.000000 2.484300 0.001400
307.000000 2.487100
Enter the values of x at which y to be calculated
301
The value of y at x=301.000000 is : 2.478597
Percent Error = 0.001236 %

Comments: Since the obtained the required value by the code is almost the same as we found
earlier by manual calculation. So this code could be implemented to solve other algebraic or
transcendental functions to find the required value using Newton’s divided difference formula.

Problem: 9 3.(a) Program to solve a system of linear equations (3 variables) with example Gauss
elimination method.

Example: use Gauss elimination to solve the system

2x + y + z = 10

3x +2 y + 3z = 18

x + 4 y + 9z = 16.
Solution: We first eliminate x from the second and third equations. For this we multiply the
first
equation by (−3/2) and add to the second to get
y + 3z = 6 (1)
Similarly, we multiply the first equation by (−1/2) and add it to the third to get
7y + 17z = 22 (2)
We thus have eliminated x from the second and third equations. Next, we have to
eliminate y from (1) and (2). For this we multiply (1) by −7 and add to (2). This gives
− 4𝑧 = −20 𝑜r 𝑧=5
The upper triangular form is therefore given by
2𝑥 + 𝑦 + 𝑧 = 1
𝑦 + 3𝑧 = 6
𝑧=6
The required solution is x = 7 , y = -9 and z = 5.

Program:

#include<iostream>
#include<iomanip>
#include<cmath>
using namespace std;
int main()

16
{
int n,i,j,k;
cout.precision(4);
cout.setf(ios::fixed);
a: cout<<"Enter the size of system"<<endl;
cin>>n;
double a[n][n+1],x[n];
cout<<"Enter the elements of system co efficients row
wise"<<endl;
for(i=0;i<n;i++)
for(j=0;j<=n;j++)
cin>>a[i][j];
//pivotisation
for(i=0;i<n;i++)
for(k=i+1;k<n;k++)
if(abs(a[i][i])<abs(a[k][i]))
{
for(j=0;j<=n;j++)
{
double temp =a[i][j];
a[i][j]=a[k][j];
a[k][j]=temp;
}
}
cout<<"The matrix after pivotaisation"<<endl;
for(i=0;i<n;i++)
{
for(j=0;j<=n;j++)
cout<<a[i][j]<<" ";
cout<<endl;
}
//Gaussian elimination
for(i=0;i<n;i++)
{
for(k=i+1;k<n;k++)
{
double t=a[k][i]/a[i][i];

for(j=0;j<=n;j++)
a[k][j]=a[k][j]-t*a[i][j];
}
}
cout<<"The new matrix after gauss elimination"<<endl;

17
 for(i=0;i<n;i++)
{
for(j=0;j<=n;j++)
cout<<a[i][j]<<setw(16);
cout<<"\n";
}
for(i=n-1;i>=0;i--)
if(a[i][i]==0)
{
cout<<"Inconsistency detected, calculation another
system"<<endl;
goto a;
}
//back substitution
for(i=n-1;i>=0;i--)
{

x[i]=a[i][n];
for(j=i+1;j<n;j++)
{
if(j!=i)
x[i]=x[i]-a[i][j]*x[j];
}
x[i]=x[i]/a[i][i];
}
cout<<"Values"<<endl;
for(i=0;i<n;i++)
cout<<x[i]<<endl;
}
Output:
Enter the size of system
3
Enter the elements of system co efficients row wise
2 1 1 10
3 2 3 18
1 4 9 16
The matrix after pivotaisation
3.0000 2.0000 3.0000 18.0000
1.0000 4.0000 9.0000 16.0000
2.0000 1.0000 1.0000 10.0000
The new matrix after gauss elimination
3.0000 2.0000 3.0000 18.0000
0.0000 3.3333 8.0000 10.0000

18
 0.0000 0.0000 -0.2000 -1.0000
Values
7.0000
-9.0000
5.0000
Comments: Since the value of the variable x,y ,z is obtained by the code is almost the same
as the value that we found earlier by manual calculation. So this code could be implemented
to solve other systems of linear equations using Gauss elimination method.

Problem: 10 3.(b) Program to solve a system of linear equations (3 variables) with example Jacobi
iterative method.

Example: Solve the system using Jacobi iterative method

10x – 5y − 2z = 3

4x − 10y + 3z = −3

x + 6y + 10z = −3
Solution: The given system can be written as
1
𝑥 = 10 (5𝑦 + 2𝑧 + 3)
1
𝑦= (4𝑥 + 3𝑧 + 3)
10
1
𝑧= (−𝑥 − 6𝑦 − 3)
10
(0) (0) (0)
Let x = 0 , y = 0 , z = 0 as initial values,
Therefore the values of x , y , z are

nth iteration x y z
1st 0.3 0.3 -0.3
2nd 0.39 0.33 -0.51
3rd 0.363 0.303 -0.537
4th 0.3441 0.2841 -0.5181
5th 0.3384 0.2822 -0.5049
6th 0.3401 0.2839 -0.5032
7th 0.3413 0.2851 -0.5044

Program:

#include<iostream>
#include<iomanip>
#include<cmath>
using namespace std;
int main()
{
cout.precision(4);
cout.setf(ios::fixed);
int k,i,j,n;

19
 cout<<"\nEnter the no. of equations\n";
cin>>n;
double a[n][n+1];
double x[n];
double x1[n];
double r[n];
double eps,temp;
cout<<"\n Enter the elements of the augmented matrix row-wise:\n";
for (i=0;i<n;i++)
for (j=0;j<=n;j++)
cin>>a[i][j];
cout<<"\n The matrix you have entered is\n";
for (i=0;i<n;i++)
{
for (j=0;j<=n;j++)
cout<<a[i][j]<< setw(10);
cout<<"\n";
}
cout<<"\n";
cout<<"\n Enter the initial values of the variables:\n";
for (i=0;i<n;i++)
cin>>x[i];
cout<<"\n Enter the accuracy upto which you want the solution:\n";
cin>>eps;
for (i=0;i<n-1;i++)
for (k=i+1;k<n;k++)
if (abs(a[i][i])<abs(a[k][i]))
for (j=0;j<=n;j++)
{
temp=a[i][j];
a[i][j]=a[k][j];
a[k][j]=temp;
}
cout<<"\n The matrix after pivoting\n";
cout<<"\n----------------------------------------------------------------------\n";
for (i=0;i<n;i++)
{
for(j=0;j<=n;j++)
{
cout<<a[i][j]<<setw(10);
}
cout<<"\n";
}
cout<<"\n";
cout<<"\n----------------------------------------------------------------------";
do
{
for (i=0;i<n;i++)
{
x1[i]=x[i];
r[i]=(a[i][n]/a[i][i]);

20
 for (j=0;j<n;j++)
{
if (j!=i)
r[i]=r[i]-(a[i][j]/a[i][i]*x[j]);
}
}
for (i=0;i<n;i++)
{
x[i]=r[i];
}
} while (fabs(x1[0]-x[0])>eps||fabs(x1[1]-x[1])>eps||fabs(x1[2]-x[2])>eps);
cout<<"\n The solution is as follows:\n";
for (i=0;i<n;i++)
cout<<"\n x["<<i<<"]="<<x[i];
return 0;
}

Output:

Enter the no. of equations

3

Enter the elements of the augmented matrix row-wise:

10 -5 -2 3
4 -10 3 -3
1 6 10 -3

The matrix you have entered is

10.0000 -5.0000 -2.0000 3.0000
4.0000 -10.0000 3.0000 -3.0000
1.0000 6.0000 10.0000 -3.0000

Enter the initial values of the variables:

0 0 0

Enter the accuracy upto which you want the solution:

.001

The matrix after pivoting

----------------------------------------------------------------------
10.0000 -5.0000 -2.0000 3.0000
4.0000 -10.0000 3.0000 -3.0000
1.0000 6.0000 10.0000 -3.0000

----------------------------------------------------------------------
The solution is as follows:

x[0] = 0.3417

21
 x[1] = 0.2852
x[2] = -0.5052

comments: Since the value of the variable x,y ,z is obtained by the code is almost the same as the
value that we found earlier by manual calculation. So this code could be implemented to solve other
systems of linear equations using the Jacobi iterative method.

Problem: 11 3.(c) Program to solve a system of linear equations (3 variables) with example
Gauss Seidal iterative method.

Example: Solve the system using Gauss –Seidal method

10x − 5y − 2z = 3
4x − 10𝑦 + 3𝑧 = −3
x + 6y + 10z = −3
Solution:
The system can be written as,
1
x = (5y + 2z + 3)
10
1
y= (4𝑥 + 3𝑧 + 3)
10
1
z= (−𝑥 − 6𝑦 − 3)
10
Taking y =0, z =0, we have ,x(1) =0.3
(0) (0)

nth iteration x y z
1st 0.3 0.42 -0.582
2nd 0.3936 0.2828 -0.5090
3rd 0.3396 0.2831 -0.5038
4th 0.3408 0.2851 -0.5052
5th 0.3416 0.2851 -0.5052
6th 0.3415 0.2850 -0.5052

Program:

#include<iostream>
#include<iomanip>
#include<cmath>
using namespace std;
int main()
{
cout.precision(7);
cout.setf(ios::fixed);
int n,i,j,k,flag=0,count=0;
cout<<"\nEnter the no. of equations\n";
cin>>n;
double a[n][n+1];
double x[n];
double eps,y;
cout<<"\nEnter the elements of the augmented matrix row-wise:\n";
for (i=0;i<n;i++)

22
 for (j=0;j<=n;j++)
cin>>a[i][j];
cout<<"\nEnter the initial values of the variables:\n";
for (i=0;i<n;i++)
cin>>x[i];
cout<<"\nEnter the accuracy upto which you want the solution:\n";
cin>>eps;
for (i=0;i<n;i++)
for (k=i+1;k<n;k++)
if (abs(a[i][i])<abs(a[k][i]))
for (j=0;j<=n;j++)
{
double temp=a[i][j];
a[i][j]=a[k][j];
a[k][j]=temp;
}
cout<<"Iter"<<setw(15);
for(i=0;i<n;i++)
cout<<"x"<<i<<setw(15);
cout<<"\n----------------------------------------------------------------------";
do
{
cout<<"\n"<<count+1<<"."<<setw(15);
for (i=0;i<n;i++)
{
y=x[i];
x[i]=a[i][n];
for (j=0;j<n;j++)
{
if (j!=i)
x[i]=x[i]-a[i][j]*x[j];
}
x[i]=x[i]/a[i][i];
if (abs(x[i]-y)<=eps)
flag++;
cout<<x[i]<<setw(15);
}
cout<<"\n";
count++;
}while(flag<n);
cout<<"\n The solution is as follows:\n";
for (i=0;i<n;i++)
cout<<"x"<<i<<" = "<<x[i]<<endl;
return 0;
}
Output:
Enter the no. of equations
3
Enter the elements of the augmented matrix row-wise:
10 -5 -2 3
4 -10 3 -3

23
 1 6 10 -3

Enter the initial values of the variables:

0 0 0

Enter the accuracy upto which you want the solution:

0.001
Iter x0 x1 x2
----------------------------------------------------------------------
1. 0.3000000 0.4200000 -0.5820000

2. 0.3936000 0.2828400 -0.5090640

3. 0.3396072 0.2831237 -0.5038349

4. 0.3407949 0.2851675 -0.5051800

5. 0.3415477 0.2850651 -0.5051938

The solution is as follows:

x0 = 0.3415477
x1 = 0.2850651
x2 = -0.5051938
Comments: Since the value of the variable x,y ,z is obtained by the code is almost the same as the
value that we found earlier by manual calculation. So this code could be implemented to solve other
systems of linear equations using the Gauss seidal iterative method.

Problem: 12 4.(a) Program to find functional value, 1st derivative, 2nd derivative and 3rd derivative
at a point from tabulated data and find the perentage error.

Example: Fit a second degree polynomial y = a + bx + cx2 to the given data (1, 0.63), (3, 2.05), (4,
4.08),
(6, 10.78) and also find functional value, 1st derivative, 2nd derivative and 3rd derivative at x = 3.
Solution:
The table of values is:

x y X2 X3 X4 xy X2y
1 0.63 1 1 1 0.63 0.63
3 2.05 9 27 81 6.15 18.45
4 4.08 16 64 256 16.32 65.28
6 10.78 36 216 1296 64.68 388.08
∑x=14 ∑y=17.54 ∑x2=62 ∑x3=308 ∑x4=1634 ∑xy=87.78 ∑x2y=472.44
Let,
y = a + bx + cx2
The normal equations are

24
 ∑y = ma + b∑x + c∑x 2

∑xy = a∑x + b∑x 2 + 𝑐∑x 3

∑x 2 𝑦 = 𝑎∑x 2 + 𝑏∑x 3 + 𝑐∑x 4

Here m = 4
Then, 4a + 14b + 62c = 17.54
14a + 62b + 308c = 87.78
62a + 308b + 1634c = 472.44
By solving , a = 1.24 , b = -1.05 , c = 0.44
The required second degree polynomial is,
y = 1.24 − 1.05x + 0.44x 2 ------------------(1)
The functional value at x=3 is,
y(3) = 1.24 − 1.05(3) + 0.44(3)2 = 2.05
Frist derivative of equation (1) is ,
𝑑𝑦 𝑑𝑦
𝑑𝑥
= 0.88𝑥 − 1.05 ; [ 𝑑𝑥 ]𝑥=3 = 0.88(3) − 1.05 = 1.59

Second derivative of equation (1) is,

𝑑2 𝑦 𝑑2 𝑦
𝑑𝑥 2
= 0.88 ; [ 𝑑𝑥 2 ]𝑥=3 = 0.88
Third derivative of equation (1) is,
𝑑3 𝑦 𝑑3 𝑦
𝑑𝑥 3
=0; [ 𝑑𝑥 3 ]𝑥=3 = 0
The percentage error is,
2.05−2.05
=| |× 100% = 0.00%
2.05

Program:

#include<iostream>
#include<iomanip>
#include<cmath>
using namespace std;
int main()
{
int n,i,j,k, N;
cout.precision(5);
cout.setf(ios::fixed);
cout<<"\nEnter the number of data points to be entered\n";
cin>>N;
double x[N],y[N];
cout<<"\nEnter the x-axis values:\n";
for (i=0;i<N;i++)
cin>>x[i];
cout<<"\nEnter the y-axis values:\n";
for (i=0;i<N;i++)

25
 cin>>y[i];
cout<<"\nEnter the degree of the polynomial you want to fit with\n";
cin>>n;
double a[n+1], B[n+1][n+2], C[n+1], X[2*n+1];
for (i=0;i<2*n+1;i++)
{
X[i]=0;
for (j=0;j<N;j++)
X[i]=X[i]+pow(x[j],i);
}
for (i=0;i<n+1;i++)
for (j=0;j<n+1;j++)
B[i][j]=X[i+j];
for(i=0;i<n+1;i++)
{
C[i]=0;
for(k=0;k<N;k++)
C[i]+=pow(x[k],i)*y[k];
}
//Now let's form the augmented matrix B by merging the column vector C
for (i=0;i<n+1;i++)
B[i][n+1]=C[i];
//Now let's set n=n+1 to make the calculation as the one used in Gaussian Elimination
n=n+1;
for (i=0;i<n;i++)
for (k=i+1;k<n;k++)
if (abs(B[i][i])<abs(B[k][i]))
for (j=0;j<=n;j++)
{
double temp=B[i][j];
B[i][j]=B[k][j];
B[k][j]=temp;
}
for (i=0;i<n-1;i++)
for (k=i+1;k<n;k++)
{
double t=B[k][i]/B[i][i];
for (j=0;j<=n;j++)
B[k][j]=B[k][j]-t*B[i][j];
}
for (i=n-1;i>=0;i--)
{
a[i]=B[i][n];
for (j=i+1;j<n;j++)
if (j!=i)
a[i]=a[i]-B[i][j]*a[j];
a[i]=a[i]/B[i][i];
}
cout<<"\nThe values of coefficients are:\n";
for (i=0;i<n;i++)
cout<<"a"<<i<<"="<<a[i]<<endl;

26
 cout<<"\nAnd the fitted Polynomial is:\ny=";
for (i=0;i<n;i++)
cout<<" + ("<<a[i]<<")"<<"x^"<<i;
cout<<"\n";
double t,sum=a[0];
cout<<"\n Functional Value and consicutive derivative";
cout<<"\n Functional Value \n";
cin>>t;
for(i=1;i<=n;i++)
{
sum=sum+a[i]*pow(t,i);
}
cout<<"FV="<<sum<<endl;
//calculation of first derivative
double sum1=a[1], prod=1;
for(i=2;i<=n;i++)
{
sum1=prod*sum1+i*a[i]*pow(t,i-1);
}
cout<<"The 1st derivative is="<<sum1<<endl;
//calculation of second derivative
double sum2 = a[2]*2;
prod = 1.0;
for (int i = 3; i <= n; i++) {
sum2 = prod * sum2 + i * (i - 1) * a[i] * pow(t, i - 2);
}
cout << "The 2nd derivative is = " << sum2 << endl;
//calculation of third derivative
double sum3 = a[3]*3*2;
prod = 1.0;
for (int i = 4; i <= n; i++) {
sum3 = prod * sum3 + i * (i - 1) * (i - 2) * a[i] * pow(t, i - 3);
}
cout << "The 3rd derivative is = " << sum3 << endl;
return 0;
}

Output:
Enter the number of data points to be entered
4
Enter the x-axis values:
1 3 4 6

Enter the y-axis values:

0.63 2.05 4.08 10.78
Enter the degree of the polynomial you want to fit with
2

27
The values of coefficients are:
a0=1.24000
a1=-1.05000
a2=0.44000
And the fitted Polynomial is:
y= + (1.24000)x^0 + (-1.05000)x^1 + (0.44000)x^2
Functional Value and consicutive derivative
Functional Value
3
FV = 2.05000
The 1st derivative is = 1.59000
The 2nd derivative is = 0.88000
The 3rd derivative is = 0.00000

Comments: The code successfully calculates the functional value, 1st derivative, 2nd derivative and 3rd
derivative at a point and matches the result obtained from manual calculations. The percentage
errors is zero.

𝐛
Problem: 13 4.(b) Numerical integration of the type∫𝐚 𝐟(𝐱) 𝐝𝐱 by using trapezoidal formula.

π/2
Example: Evaluate ∫0 sin x dx using Trapezoidal formula.
Solution:
Let h = π/12= 3.1416/12
x y= sinx
x0 = 0 y0= 0
x1= π/2 y1= 0.25882
x2= 2π/12 y2= 0.5
x3= 3 π/12 y3= 0.70711
x4= 4 π/12 y4= 0.86603
x5= 5 π/12 y5= 0.96593
X6= 6 π/12 y6= 1.0

By Trapezoidal rule we have ,

π
h
∫02 sin x dx = 2 [(y0 + y6) + 2(y1 + y2 + y3 + y4 + y5)]

π
= [0 + 1 + 2(0.25882 + 0.5 + 0.70711 + 0.86603 + 0.96593)]
24

= 0.9943

Program:
#include<iostream>

28
 #include<cmath>
using namespace std;
double f(double x)
{
double a=sin(x);
return a;
}
int main()
{
int n,i;
double a,b,h,sum=0,integral;
cout<<"Enter the limits of integration,\nInitial limit,a=";
cin>>a;
cout<<"Final limit, b=";
cin>>b;
cout<<"Enter the no. of subintervals, n=";
cin>>n;
double x[n+1],y[n+1];
h=(b-a)/n;
for (i=0;i<=n;i++)
{
x[i]=a+i*h;
y[i]=f(x[i]);
}
for (i=1;i<n;i++)
{
sum=sum+y[i];
}
integral=h/2.0*((y[0]+y[n])+2.0*sum);
cout<<"The definite integral is "<<integral<<endl;
return 0;
}

Output:
Enter the limits of integration,
Initial limit,a =0
Final limit, b =1.5708
Enter the no. of subintervals, n =6
The definite integral is 0.994286

Comments: The code successfully calculates the value of an integration using trapezoidal formula
and matches the result obtained from manual calculations. This indicates that the code can be
applied to solve other integration problems as well. The key takeaway is that the trapezoidal formula
implementation in the code appears to be accurate and can be reused for a variety of integration
tasks.

𝐛
Problem: 14 4.(c) Numerical integration of the type∫𝐚 𝐟(𝐱) 𝐝𝐱 by using Simson’s 1/3 formula.
π/2
Example: Evaluate ∫0 sin x dx using Simson’s 1/3 formula.
Solution:
Let h = π/12 =3.1416/12

29
 x y= sinx
x0 = 0 y0= 0
x1= π/2 y1= 0.25882
x2= 2π/12 y2= 0.5
x3= 3 π/12 y3= 0.70711
x4= 4 π/12 y4= 0.86603
x5= 5 π/12 y5= 0.96593
X6= 6 π/12 y6= 1.0

By Simson’s 1/3 rule we have ,

π
h
∫02 sin x dx = 3 [(y0 + y6) + 4(y1 + y3 + y5) + 2(𝑦2 + 𝑦4)]

π
= [0 + 1 + 4(0.25882 + 0.70711 + 0.96593) + 2(0.5 + 0.86603)]
36

= 1.000032

Program:

#include<iostream>
#include<cmath>
using namespace std;
double f(double x)
{
double a=sin(x);
return a;
}
int main()
{
cout.precision(5);
cout.setf(ios::fixed);
int n,i;
double a,b,h,sum=0,integral;
cout<<"\nEnter the interval of integration,\n\nLower limit a= ";
cin>>a;
cout<<"\nUpper limit, b=";
cin>>b;
cout<<"\nEnter the no. of subintervals\n\n";
cout<<"\nPlease enter a number that is multiple of 2:";
cin>>n;
double x[n+1],y[n+1];
h=(b-a)/n;
for (i=0;i<=n;i++)
{
x[i]=a+i*h;
y[i]=f(x[i]);
}
for (i=1;i<n;i+=2)
{
sum=sum+4.0*y[i];

30
 }
for (i=2;i<n;i+=2)
{
sum=sum+2.0*y[i];
}
integral=h/3.0*(y[0]+y[n]+sum);
cout<<"\nThe definite integral is "<<integral<<"\n"<<endl;
return 0;
}

Output:
Enter the interval of integration,
Lower limit a= 0
Upper limit, b=1.5708
Enter the no. of subintervals
Please enter a number that is multiple of 2: 6
The definite integral is 1.00003

Comments: The code successfully calculates the value of an integration using Simpson's 1/3 rule
and matches the result obtained from manual calculations. This indicates that the code can be
applied to solve other integration problems as well. The key takeaway is that the Simpson's 1/3 rule
implementation in the code appears to be accurate and can be reused for a variety of integration
tasks.

𝐛
Problem: 15 4.(d) Numerical integration of the type∫𝐚 𝐟(𝐱) 𝐝𝐱 by using Simson’s 3/8 formula.
π/2
Example: Evaluate ∫0 sin x dx using Simson’s 3/8 formula.
Solution:
Let h = π/12 =3.1416/12
x y= sinx
x0 = 0 y0= 0
x1= π/2 y1= 0.25882
x2= 2π/12 y2= 0.5
x3= 3 π/12 y3= 0.70711
x4= 4 π/12 y4= 0.86603
x5= 5 π/12 y5= 0.96593
X6= 6 π/12 y6= 1.0

By Simson’s 3/8 rule we have ,

π
3h
∫02 sin x dx = [(y0 + y6) + 3(y1 + y2 + y4 + y5) + 2y3]
8

3π
= [0 + 1 + 3(0.25882 + 0.5 + 0.86603 + 0.96593) + 2(0.70711)]
96

= 1.000065
Program:

#include<iostream>
#include<cmath>
using namespace std;

31
 #define f(x) (sin(x))
int main()
{
cout.precision(6);
cout.setf(ios::fixed);
int n,i;
double a,b,c,h,sum=0,integral;
cout<<"\nEnter the interval of integration,\n\nLower limit,a= ";
cin>>a;
cout<<"\nUpper limit, b="; //get the limits of integration
cin>>b;
cout<<"\nEnter the no. of subintervals\n\n";
cout<<"\nPlease enter a number that is multiple of 3:\n\n";
cin>>n;
double x[n+1],y[n+1];
h=(b-a)/n;
for (i=0;i<=n;i++){
x[i]=a+i*h;
y[i]=f(x[i]);
}
for (i=1;i<n;i++)
{
if (i%3==0)
sum=sum+2*y[i];
else
sum=sum+3*y[i];
}
integral=3*h/8*(y[0]+y[n]+sum);
cout<<"\nThe definite integral is "<<integral<<"\n"<<endl;
return 0;
}

Output:
Enter the interval of integration,
Lower limit,a= 0
Upper limit, b= 1.5708
Enter the no. of subintervals
Please enter a number that is multiple of 3:
6

The definite integral is 1.000063

Comments : The code successfully calculates the value of an integration using Simpson's 3/8 rule
and matches the result obtained from manual calculations. This indicates that the code can be
applied to solve other integration problems as well. The key takeaway is that the Simpson's 3/8 rule
implementation in the code appears to be accurate and can be reused for a variety of integration
tasks.

32
 𝐝𝐱
Problem: 16 5.(a) Numerical solution of first order DE of the type = 𝐟(𝐭, 𝐱) by Euler formula at
𝐝𝐭
a point and for particular range with example.

Example: Solve the following DE using Euler method

dy
dx
= − y with the condition y(0) = 1
Find the value of y at x= 0.04 .
Solution: we know the general formula is
yn+1 = yn + hf(x1 , y1 ) where n = 1 ,2 ,3 ,………
Let h = 0.01
y(0.01) = 1 + 0.1(−1) = 0.99
y(0.02) = 0.99 + 0.01(−0.99) = 0.9801
y(0.03) = 0.9801 + 0.01(−0.9801) = 0.9703
y(0.04) = 0.9703 + 0.01(−0.9703) = 0.9606
The exact solution is y =𝑒 −𝑥 and from this the value at x = 0.04 is 0.9608
Program:

#include<iostream>
#include<iomanip>
#include<cmath>
using namespace std;
#define f(x,y) (-y)
int main()
{
int n,i;
double x0,y0,h,fx;
cout.precision(5);
cout.setf(ios::fixed);
cout<<"\nEnter the initial condition\n";
cout<<"\nEnter the initial value of x\n";
cin>>x0;
cout<<"\nEnter the initial value of y\n";
cin>>y0;
cout<<"\nFor what value of x do you want to find the value of y\n";
cin>>fx;
cout<<"\nEnter the number of subdivisions\n";
cout<<"\nPlease enter a large value of \n";
cin>>n;
h=(fx-x0)/n;
double x[n],y[n];
x[0]=x0;
y[0]=y0;
cout<<" x"<<setw(10)<<"y(x)";
cout<<"\n-----------------------\n";
for(i=0;i<n;i++)
{
y[i+1]=y[i]+h*f(x[i],y[i]);
x[i+1]= x[i]+h;
cout<<x[i+1]<<setw(10)<<y[i+1]<<endl;
}

cout<<"y("<<x[n]<<") = "<<y[n]<<endl;

33
 double per_err;
per_err=(abs((exp(-fx))-y[n]))/(exp(-fx))*100;
cout<<"Percentage error = "<<per_err<<" % "<<endl;
return 0;
}
Output:
Enter the initial condition
Enter the initial value of x
0
Enter the initial value of y
1
For what value of x do you want to find the value of y
0.04
Enter the number of subdivisions
Please enter a large value of
5
x y(x)
--------------------------------
0.00800 0.99200
0.01600 0.98406
0.02400 0.97619
0.03200 0.96838
0.04000 0.96063
y(0.04000) = 0.96063
Percentage error = 0.01608 %

Comments: Since the value of y at x = 0.04 obtained by the code is same as the value that
we found earlier by manual calculation. And the Percent Error is just 0.01608 % . It is very
small .So this code could be implemented to solve other 1st order differential equations by
Euler's method.

𝐝𝐱
Problem: 17 5.(b) Numerical solution of first order DE of the type = 𝐟(𝐭, 𝐱) by Modified Euler
𝐝𝐭
formula at a point and for particular range with example.

Example: Let us consider the first order DE with IVP,

dy
= y + x4 y(0) = −23
dx
Using Modified Euler’s method find the value of y at x =1 .
Solution:
Here f(x,y)= y + x4
Let h = 0.2
Modified Euler’s formula,
h
y1 = y0 + [f(x0 , y0 ) + f(x1 , y1 )]
2
Therefore,
0.2
y1 (0.2) = −23 + [𝑓(x0 , y0 ) + f(x1 , y1 )]
2

= −23 + 0.1[y0 + 𝑥04 + y0 + 0.2(y0 + 𝑥04 )] = −28.05984

Similarly,
y2 (0.4) = −33.66924

34
 𝑦3 (0.6) = −40.3899776
𝑦4 (0.8) = −48.4244531
𝑦5 (1.0) = −57.9963837
The value at x =1 is -57.9963837

Program:
#include<iostream>
#include<iomanip>
#include<cmath>
using namespace std;
#define f(x,y) (y+x*x*x*x)
int main()
{
int n,i;
double x0,y0,h,X;
cout.precision(7);
cout.setf(ios::fixed);
cout<<"\nEnter the initial condition\n";
cout<<"\nEnter the initial value of x\n";
cin>>x0;
cout<<"\nEnter the initial value of y\n";
cin>>y0;
cout<<"\nFor what value of x do you want to find the value of y\n";
cin>>X;
cout<<"\nEnter the number of subdivisions\n";
cout<<"\nPlease enter a large value of \n";
cin>>n;
h=(X-x0)/n;
double x[n],y1[n],y[n];
x[0]=x0;
y[0]=y0;
y1[0]=0;
cout<<" x"<<setw(10)<<"y(x)";
cout<<"\n-----------------------\n";
for(i=0;i<n;i++)
{
y[i+1]=y[i]+h*f(x[i],y[i]);
x[i+1]=x[i]+h;
y1[i+1]=y[i]+h/2.0*(f(x[i],y[i])+f(x[i+1],y[i+1]));
cout<<x[i+1]<<setw(10)<<y1[i+1]<<endl;
}
cout<<"y("<<x[n]<<") = "<<y1[n]<<endl;
return 0;
}

Output:
Enter the initial condition
Enter the initial value of x
0

35
Enter the initial value of y
-23
For what value of x do you want to find the value of y
1
Enter the number of subdivisions
Please enter a large value of
5
x y(x)
---------------------------------------
0.2000000 -28.0598400
0.4000000 -33.6692480
0.6000000 -40.3899776
0.8000000 -48.4244531
1.0000000 -57.9963837
y(1.0000000) = -57.9963837
Comments: Since the value y at x = 1.0 obtained by the code is same as the value that we found
earlier by manual calculation.So this code could be implemented to solve other 1st order differential
equations by Euler's Modified method.
𝐝𝐱
Problem: 18 5.(c) Numerical solution of first order DE of the type 𝐝𝐭 = 𝐟(𝐭, 𝐱) by Runge-Kutta 4th
order formula at a point and for particular range with example.
dy
Example: Given = y − x where y(0)=2 , find y(0.1) correct to four decimal places by using Runge-
dx
Kutta 4th order method.
Solution: Runge-Kutta 4th order formula is,
k1 = hf(x0 , y0 )
1 1
k 2 = hf(x0 + 2 h, y0 + 2 k1 )
1 1
k 3 = hf(x0 + 2 h, y0 + 2 k 2 )
k 4 = hf(x0 + h, y0 + k 3 )
1
y1 = y0 + 6 (k1 + 2k 2 + 2k 3 + k 4 )
To determine y(0.1),we have x0 = 0 , y0 = 2 and h= 0.1 .then we get,
K1 = 0.2
k 2 = 0.205
K 3 = 0.20525
K 4 = 0.21053
1
y(0.1) = 2 + (0.2 + 2(0.205) + 2(0.20525) + 0.21053) = 2.2052
6

Program:
#include<iostream>

36
#include<iomanip>
#include<cmath>
using namespace std;
#define f(x,y) (y-x)
int main()
{
int n,i;
double x0,y0,h,fx,k1,k2,k3,k4;
cout.precision(7);
cout.setf(ios::fixed);
cout<<"\nEnter the initial condition\n";
cout<<"\nEnter the initial value of x\n";
cin>>x0;
cout<<"\nEnter the initial value of y\n";
cin>>y0;
cout<<"\nFor what value of x do you want to find the value of y\n";
cin>>fx;
cout<<"\nEnter the number of subdivisions\n";
cin>>n;
h=(fx-x0)/n;
double x[n],y[n];
x[0]=x0;
y[0]=y0;
cout<<" x"<<setw(10)<<"y(x)";
cout<<"\n-----------------------\n";
for(i=0;i<n;i++)
{
k1=h*f(x[i],y[i]);
k2=h*f(x[i]+h/2.0, y[i]+k1/2.0);
k3=h*f(x[i]+h/2.0, y[i]+k2/2.0);
k4=h*f(x[i]+h, y[i]+k3);
y[i+1]=y[i]+(1/6.0)*(k1+2*k2+2*k3+k4);
x[i+1]=x[i]+h;
cout<<x[i+1]<<setw(10)<<y[i+1]<<endl;
}
cout<<"\n\nThe approximate value of y at x="<<x[n]<<" is="<<y[n]<<endl;
return 0;
}
Output:
Enter the initial condition
Enter the initial value of x
0
Enter the initial value of y
2
For what value of x do you want to find the value of y
0.1
Enter the number of subdivisions

37
 10
x y(x)
-----------------------------------
0.0100000 2.0201505
0.0200000 2.0404030
0.0300000 2.0607586
0.0400000 2.0812182
0.0500000 2.1017830
0.0600000 2.1224539
0.0700000 2.1432321
0.0800000 2.1641186
0.0900000 2.1851145
0.1000000 2.2062209
The approximate value of y at x=0.1000000 is=2.2062209
Comments: Since the value y at x = 0.1 obtained by the code is same as the value that we found
earlier by manual calculation.So this code could be implement to solve other 1st order differential
equation by Range-kutta 4th order method.
Problem: 19 5.(d) Numerical solution of system of 1st order DE by Runge-Kutta 4th order formula
at a point and for particular range with example.
dy dz
Example: Solve the system using Runge-Kutta 4th order method = yz + x , = xz + y , where
dx dx
y(0)=1 and z(0)= -1 .Find y(0.1) and z(0.1).
Solution:
dy
Supppose = yz + x = f(x, y, z)
dx
dz
dx
= xz + y = g(x, y, z)
we know Runge-Kutta 4th order formula is ,
Here h=0.1,
k1 = hf(x0 , y0 , z0 )
= 0.1 f(0,1, −1) = −0.1
l1 = hg(x0 , y0 , z0 )
= 0.1g(0,1, −1) = 0.1
1 1 1
k 2 = hf(x0 + 2 h, y0 + 2 k1 , z0 + 2 l1 )
= 0.1f(0.05,0.95, −0.95) = −0.085
1 1 1
l2 = hg(x0 + 2 h, y0 + 2 k1 , z0 + 2 l1 )
= 0.1g(0.05,0.95,0.95) = 0.0902
1 1 1
k 3 = hf(x0 + 2 h, y0 + 2 k 2 , z0 + 2 l2 )

38
 = −0.08641
1 1 1
l3 = hg(x0 + 2 h, y0 + 2 k 2 , z0 + 2 l2 )

= 0.09096
k 4 = hf(x0 + h, y0 + k 3 , z0 + l3 )
= −0.073048
l4 = hg(x0 + h, y0 + k 3 , z0 + l3 )
= 0.08226
1
y1 = y0 + (k1 + 2k 2 + 2k 3 + k 4 )
6
1
= 1 + 6 [−0.1 + 2(−0.085) + 2(−0.08641) − 0.073048]
= 0.9141
1
z1 = z0 + (l1 + 2l2 + 2l3 + l4 )
6
1
= −1 + 6 [0.1 + 2(0.0902) + 2(0.09096) + 0.08226]
= −0.9092
Program:

#include<iostream>
#include<cmath>
#include<iomanip>
using namespace std;
#define dy_dx(x,y,z) (z)
#define dz_dx(x,y,z) (-y)
int main()
{
cout.precision(6);
cout.setf(ios::fixed);
int i,n;
double x0,y0,z0,h,fx;
cout<<"Enter the initial value of x:\n";
cin>>x0;
cout<<"Enter the initial value of y corresponding to x:\n";
cin>>y0;
cout<<"Enter the initial value of z corresponding to x:\n";
cin>>z0;
cout<<"Enter the value of x up to which you want to find the value of y & z:\n";
cin>>fx;
cout<<"Enter the value of h\n";
cin>>h;
n=(fx-x0)/h;
double x[n+1],y[n+1],z[n+1],k1,k2,k3,k4,l1,l2,l3,l4;
x[0]=x0,
y[0]=y0,
z[0]=z0;
for(i=0;i<n;i++)
{

39
k1=h*dy_dx(x[i],y[i],z[i]);
l1=h*dz_dx(x[i],y[i],z[i]);
k2=h*dy_dx((x[i]+h/2.0),(y[i]+k1/2.0),(z[i]+l1/2.0));
l2=h*dz_dx((x[i]+h/2.0),(y[i]+k1/2.0),(z[i]+l1/2.0));
k3=h*dy_dx((x[i]+h/2.0),(y[i]+k2/2.0),(z[i]+l2/2.0));
l3=h*dz_dx((x[i]+h/2.0),(y[i]+k2/2.0),(z[i]+l2/2.0));
k4=h*dy_dx((x[i]+h),(y[i]+k3),(z[i]+l3));
l4=h*dz_dx((x[i]+h),(y[i]+k3),(z[i]+l3));
y[i+1]=y[i]+(1/6.0)*(k1+2*k2+2*k3+k4);
z[i+1]=z[i]+(1/6.0)*(l1+2*l2+2*l3+l4);
x[i+1]=x[i]+h;
}
for(i=0; i<=n; i++)
{
cout<<"y("<<x[i]<<"): "<<y[i]<<endl;
cout<<"z("<<x[i]<<"): "<<z[i]<<endl;
cout<<"\n";
}
return 0;
}

Output:
Enter the initial value of x:
0
Enter the initial value of y corresponding to x:
1
Enter the initial value of z corresponding to x:
-1
Enter the value of x up to which you want to find the value of y & z:
0.1
Enter the value of h
0.1
y(0.000000): 1.000000
z(0.000000): -1.000000

y(0.100000): 0.913936
z(0.100000): -0.909218
Comments: Since the value of y and z at x = 0.1 obtained by the code is approximately same as the
value that we found earlier by manual calculation.So this code could be implement to solve the
system of 1st order differential equation by Range-kutta 4th order method.

40
 𝛛𝐮 𝛛𝟐 𝐮
Problem: 20 6.(a) Numerical solution of parabolic type PDE 𝐝𝐭
= 𝐜𝟐 𝐝𝐱 𝟐
by finite difference
formula at a point and for particular range with example.

∂u ∂2 u
Example: Solve the heat conduction problem = 2 subject to the condition u(x, 0) = sinπx ,
dt dx
0 ≤ x ≤ 1, and u(0, t) = 0 = u(1, t). Use Bender-Schmidt’s method find the value of u(x,0.04).

Solution:
Let h=0.2 then l=λh2=0.02
Initial values of u are
u00 = 0, u10 = 0.5878, u02 = 0.9510,
u03 = 0.9510, u04 = 0.5878, u05 = 0
The Bender-Schmidt formula gives,
1 1
u11 = (0.9510) = 0.4755, u12 = (0.5878 + 0.9510) = 0.7694,
2 2
u13 = 0.7694, u14 = 0.4755
Also,

1
𝑢12 = (0.7694) = 0.3847, 𝑢22 = 0.62245, 𝑢32 = 0.62245, 𝑢42 = 0.3847
2

Program:

#include <iostream>
#include<iomanip>
#include<cmath>
using namespace std;

//Function of initial condition

double f(double x)
{
double pi=3.14159265;
double a=sin(pi*x);
return a;
}

int main()
{
cout.precision(5);
cout.setf(ios::fixed);
int i, j, n, m;
double xl, yl, xr, yr, lam, td, xd, t=0.0, T, X, pi=3.14159265;
cout<<"\nPlease enter the boundary conditions\n";
cout<<"\nThe left boundary condition is\n";
cin>>xl>>yl;
cout<<"\nThe right boundary condition is\n";
cin>>xr>>yr;
cout<<"\nEnter the number of subdivision x direction\n";
cin>>n;
a:cout<<"\nEnter the number of subdivision t direction\n";

41
cin>>m;
xd=(xr-xl)/n;
cout<<"\nEnter the final time of simulation\n";
cin>>T;
td=(T-t)/m;
lam=td/(xd*xd);
if(lam>0.5)
{
cout<<"\nPlease enter a larger number of subdivision in t direction\n";
goto a;
}
else
cout<<"\nThe given subdivisions satsify the convergence criterion of the explicit finite
difference scheme\n";
double x[n+1], tt[m+1], y[n+1][m+1];
x[0]=xl; tt[0]=0.0; x[n]=xr; tt[n]=T; y[0][0]=yl; y[n][0]=yr;
//Initialize the initial conditions at all x nodes at time t=0
for(i=1;i<n;i++)
{
x[i]=i*xd;
y[i][0]=f(x[i]);
}
//for(i=0;i<=n;i++)
// cout<<"\n"<<x[i]<<setw(20);
//Main loop for iteration
for(j=0;j<=m;j++)
{
for(i=1;i<n;i++)
{
y[i][j+1]=y[i][j]+lam*(y[i-1][j]-2*y[i][j]+y[i+1][j]);

}
y[0][j+1]=0; //effect of zero boundary conditions at all time steps
y[n][j+1]=0;
}
b: cout<<"\nPlease enter the value of t at which you would like to find of x\n";
cin>>T;
int k,l;
k=T/td+1;
cout<<"\nThe value of y at t="<<T<<" are\n\n";
for(i=0;i<=n;i++)
cout<<y[i][k]<<setw(10);
int F;
cout<<"\n\n\nDo you want to see the x values at any other value of t??\n";
cout<<"\n If yes please enter 1; if no please enter 0\n";
cin>>F;
if(F<=0)
goto c;
else
goto b;

42
 c: cout<<"\nDo you want to see the values at a particular x for all t??\n";
cout<<"\nIf yes please enter 1; if no please enter 0\n";
cin>>F;
if(F<=0)
goto d;
else
goto e;

e: cout<<"\nThe value of x at which you would like to see the x value for all t\n";
cin>>X;
//int k;
k=X/xd+1;
cout<<"\nThe value of y at x="<<X<<" are\n\n";
for(j=0;j<=m;j++)
cout<<y[k][j]<<setw(10);
cout<<"\n\n\nDo you want to see the y values at any other x for all t??\n";
cout<<"\nIf yes please enter 1; if no please enter 0\n";
cin>>F;
if(F<=0)
goto d;
else
goto e;
d:cout<<"\n\n\nDo you want to see the y value at a particular x and t??\n";
cout<<"\nIf yes please enter 1; if no please enter 0\n";
cin>>F;
if(F<=0)
cout<<"\nPlease re-run the program in case you change your mind and want to find any
other information\n";
else
{
double p,q;
cout<<"\nPlease enter the value of x and t at which you would like to see the y
value\n";
cin>>p>>q;
k=p/xd+1; l=q/td+1;
cout<<"\nThe value of y at x="<<p<<" and t="<<q<<" is\n\n";
cout<<y[k][l];
}
return 0;
}

Output:
Please enter the boundary conditions
The left boundary condition is
0
0
The right boundary condition is
1

43
 0
Enter the number of subdivision x direction
5
Enter the number of subdivision t direction
10
Enter the final time of simulation
0.04
The given subdivisions satsify the convergence criterion of the explicit finite difference scheme
Please enter the value of t at which you would like to find of x
0.04
The value of x at t=0.04000 are
0.00000 0.38297 0.61966 0.61966 0.38297 0.00000
Do you want to see the x values at any other value of t??
If yes please enter 1; if no please enter 0
0
Comments: The code successfully calculates the value using least square method and almost
matches the result obtained from manual calculations. This indicates that the code can be applied to
solve other problems as well.

𝛛𝟐 𝐮 𝛛𝟐 𝐮
Problem: 21 6.(b) Numerical solution of hyperbolic type PDE 𝐝𝐭 𝟐
= 𝐜𝟐 𝐝𝐱 𝟐
by finite difference
formula at a point and for particular range with example.

∂2 u ∂2 u
Example: Solve the Hyperabolic PDE dt2 = 4 dx2 , subject to the condition
u(0, t) = 0 = u(4, t), ut (x, 0) = 0, u(x, 0) = 4x − x 2 .
Solution:
∂2 u ∂2 u
Given that, dt2 = 4 dx2 ,
u(0, t) = 0 = u(4, t), ut (x, 0) = 0, u(x, 0) = 4x − x 2

Let , h=1 and α=1, so l=0.5

We have,
uk0 = uk4 = 0 for all k
ut (x, 0) = 0, we obtain u−1 1
i = ui
Further,
u(x, 0) = 4x − x 2
Implies that u0i = 4i − i2 , since h = 1
Then,
𝑢00 = 0
𝑢10 = 3
𝑢20 = 4
𝑢30 = 3
𝑢40 = 0

44
Program:

#include <iostream>
#include <cmath>
#include <iomanip>
using namespace std;
#define u_x_0(x) (4*x-x*x) //f(x)=(4*x-x*x) assigned as 1st initial condition
#define du_dt(x) (0) //g(x)=0 assigned as 2nd initial condition
int main()
{
cout.precision(6);
cout.setf(ios::fixed);
int i,j,m,n,option;
double c,alpha,h,k,x,t,x_i,x_f,t_i,t_f,u_i,u_f;
cout<<"Numerical solution of wave equation using finite difference formula"<<endl;
cout<<"Enter value of c"<<endl;
cin>>c;
cout<<"Enter initial x and corresponding u (1st boundary condition)"<<endl;
cin>>x_i>>u_i; //1st boundary condition u(x_i,t)=u_i
cout<<"Enter final x and corresponding u (2nd boundary condition)"<<endl;
cin>>x_f>>u_f; //2nd boundary condition u(x_f,t)=u_f
t_i=0; //initial simulation time assigned
cout<<"Enter final simulation time"<<endl;
cin>>t_f;
cout<<"Enter value of h and k"<<endl;
z:
cin>>h>>k;
alpha = c*(k/h); //calculated alpha
//check alpha limitation
if(alpha>1)
{
cout<<"please enter new value of h & k such that alpha<1 where, alpha=c*(k/h)"<<endl;
goto z;
}
//calculated division of x & t
m=(x_f-x_i)/h;
n=(t_f-t_i)/k;
//u defined as [(n+1) by (m+1)] dimensional array
double u[n+1][m+1];
//calculated 1st & last column
for(i=0; i<=n; i++)
{
u[i][0]=u_i; //1st column
u[i][m]=u_f; //last column
}
//calculated 1st row
for(j=1; j<=m-1; j++)
{
x=x_i+j*h; //transforming j into x
u[0][j]= u_x_0(x);
}

45
//calculated 2nd row
for(j=1; j<=m-1; j++)
{
x=x_i+j*h; //transforming j into x
u[1][j]=((pow(alpha,2))*(u[0][j+1]+u[0][j-1]))/2.0+(1-pow(alpha,2))*u[0][j]+k*du_dt(x);
}
//row wise calculated 3rd to last row
for(i=1; i<=n-1; i++)
{
for(j=1; j<=m-1; j++)
{
x=x_i+j*h; //transforming j into xu[i+1][j]=(pow(alpha,2))*(u[i][j1]+u[i][j+1])+2*(1-
pow(alpha,2))*u[i][j]-u[i-1][j];
}
}
//printing value of u(x,t)
while(1)
{
cout<<"\nChoose option how you want result:"<<endl;
cout<<"1. x, t both vary"<<endl;
cout<<"2. x vary, t constant"<<endl;
cout<<"3. x constant, t vary"<<endl;
cout<<"4. x, t both constant"<<endl;
cin>>option;
if(option==1)
{
for( i=0; i<=n; i++)
{
t=t_i+i*k; //transforming i into t
for( j=0; j<=m; j++)
{
x=x_i+j*h; //transforming j into xcout<<"u("<<x<<","<<t<<")= "<<u[i][j]<<endl;
}
}
}
else if(option==2)
{
cout<<"Enter value of t where you want to find the value of u:"<<endl;
cin>>t;
i=(t-t_i)/k; //transforming t into i
for( j=0; j<=m; j++)
{
x=x_i+j*h; //transforming j into x
cout<<"u("<<x<<","<<t<<")= "<<u[i][j]<<endl;
}
}
else if(option==3)
{
cout<<"Enter value of x where you want to find the value of u:"<<endl;
cin>>x;
j=(x-x_i)/h; //transforming x into j

46
 for( i=0; i<=n; i++)
{
t=t_i+i*k; //transforming i into t
cout<<"u("<<x<<","<<t<<")= "<<u[i][j]<<endl;
}
}
else if(option==4)
{
cout<<"Enter value of x where you want to find the value of u:"<<endl;
cin>>x;
cout<<"Enter value of t where you want to find the value of u:"<<endl;
cin>>t;
i=(t-t_i)/k; //transforming t into i
j=(x-x_i)/h; //transforming x into j
cout<<"u("<<x<<","<<t<<")= "<<u[i][j]<<endl;
}
}
return 0;
}

Output:
Numerical solution of wave equation using finite difference formula
Enter value of c
2
Enter initial x and corresponding u (1st boundary condition)
0
0
Enter final x and corresponding u (2nd boundary condition)
4
0
Enter final simulation time
1
Enter value of h and k
1
0.5

Choose option how you want result:

1. x, t both vary
2. x vary, t constant
3. x constant, t vary
4. x, t both constant
2
Enter value of t where you want to find the value of u:
0
u(0.000000,0.000000)= 0.000000
u(1.000000,0.000000)= 3.000000
u(2.000000,0.000000)= 4.000000
u(3.000000,0.000000)= 3.000000
u(4.000000,0.000000)= 0.000000

Comments: The code successfully calculates the value and matches the result obtained from manual

47
calculations. This indicates that the code can be applied to solve other problems as well.

𝛛𝟐 𝐮 𝛛𝟐 𝐮
Problem: 22 6.(c) Numerical solution of elliptic type PDE 𝐝𝐱 𝟐
+ 𝐝𝐲𝟐 = 𝟎 by finite difference
formula at a point and for particular range with example.

∂2 u ∂2 u
Example: Solve the equation + 2 = 0 in the domain of figure below
dx2 dy

1 1

0 0
𝑢4 𝑢3

0 0
𝑢1 𝑢2

0 0
We obtain the approximate values of u1, u2, u3, u4 as follows,
(1) 1
𝑢1 = (0 + 0 + 0 + 1) = 0.25
4
(1) 1
𝑢2 = 4 (0 + 0 + 0 + 1) = 0.25
(3) 1
𝑢1 = 4 (1 + 1 + 0 + 0) = 0.5
(4) 1
𝑢1 = 4 (1 + 1 + 0 + 0) = 0.5
The iterations have been continued using Jacobi method,and seven successive iterates are given
below:

U1 U2 U3 U4
0.1875 0.1875 0.4375 0.4375
0.15625 0.15625 0.40625 0.40625
0.14062 0.14062 0.39062 0.39062
0.13281 0.13281 0.38281 0.38281
0.12891 0.12891 0.37891 0.37891
0.12695 0.12695 0.37695 0.37695
0.12598 0.12598 0.37598 0.37598

Program:

#include <iostream>
#include <cmath>
#include <vector>
using namespace std;
double tolerance = 1e-6; // Convergence tolerance
int maxIterations = 10000; // Maximum number of iterations

48
void jacobiMethod(vector<double>& u) {
vector<double> u_old(4, 0.0); // Old values for comparison
int iteration = 0;
while (iteration < maxIterations) {
u_old = u;
// Jacobi update
u[0] = 0.25 * (0 + u_old[1] + 0 + u_old[3]);
u[1] = 0.25 * (u_old[0] + 1 + 0 + u_old[2]);
u[2] = 0.25 * (u_old[1] + 1 + 0 + u_old[3]);
u[3] = 0.25 * (u_old[0] + 0 + u_old[2] + 0);
// Check for convergence
double maxDiff = 0.0;
for (int i = 0; i < 4; ++i) {
maxDiff = max(maxDiff, fabs(u[i] - u_old[i]));
}
if (maxDiff < tolerance) {
break;
}
iteration++;
}
}
int main() {
vector<double> u(4, 0.0); // Initial guess for u1, u2, u3, u4
jacobiMethod(u);
cout<< "The values of u are \n ";
cout << "u1: " << u[0] << endl;
cout << "u2: " << u[3] << endl;
cout << "u3: " << u[2] << endl;
cout << "u4: " << u[1] << endl;
return 0;
}

Output:
The values of u are
u1: 0.124999
u2: 0.124999
u3: 0.374999
u4: 0.374999

Comments: The code almost successfully calculates the value and near the result obtained from
manual calculations. But this program is not generalized for other problems.

Problem: 23 7.(a) Curve fitting of the form y = a + bx .

Example: The values of x and y are given , (0,-1), (2,5), (5,12), (7,20) fit a curve to the above data.
Solution:
The table of values is given below:

x y X2 xy
0 -1 0 0
2 5 4 10

49
 5 12 25 60
7 20 49 140
∑x=14 ∑y=36 ∑x2=78 ∑xy=210

The normal equations are,

∑y = ma + b∑x
∑xy = a∑x + b∑x 2
Here m = 4
And
4a0 + 14a1 = 36
14a0 + 78a1 = 210
Solving these two equations we get,
a0 = −1.1381 and a1 = 2.8966
Hence the best straight line fit is given by,
y = −1.1381 + x(2.8966)

Program:

#include <iostream>
using namespace std;
int main()
{
int n; // number of data points
cout << "Enter the number of data points: ";
cin >> n;
double **data = new double*[n];
for (int i = 0; i < n; i++) {
data[i] = new double[2]; // each data point has x and y coordinates
}

// Input the data points from the user

cout << "Enter the data points (x y):" << endl;
for (int i = 0; i < n; i++)
{
cout << "Data point " << i+1 << ": ";
cin >> data[i][0] >> data[i][1]; // x-coordinate and y-coordinate
}
// Compute the sums needed for calculating the slope and intercept of the line
double sumX = 0, sumY = 0, sumXY = 0, sumX2 = 0;
for (int i = 0; i < n; i++) {
sumX += data[i][0];
sumY += data[i][1];
sumXY += data[i][0] * data[i][1];
sumX2 += data[i][0] * data[i][0];
}
double slope = (n * sumXY - sumX * sumY) / (n * sumX2 - sumX * sumX);
double intercept = (sumY - slope * sumX) / n;
// Output the equation of the fitted line
cout << "The equation of the fitted line is: y = " << intercept<<"+ " << slope << "x " endl;
return 0;
}

50
Output:
Enter the number of data points: 4
Enter the data points (x y):
Data point 1: 0 -1
Data point 2: 2 5
Data point 3: 5 12
Data point 4: 7 20
The equation of the fitted line is : y = -1.13793 + 2.89655 x
Comments: Since the fitted polynomial obtained by the code is same as the fitted polynomial that
we found earlier by manual calculation. So this code could be implement to fitted other polynomial
of degree one for a given set of points.

Problem: 24 7.(b) Curve fitting of the form y = a + bx + cx2.

Example: Fit a second degree parabola y = a + bx + cx2 to the given data (1, 0.63), (3, 2.05), (4, 4.08),
(6, 10.78).
Solution:
The table of values is:

x y X2 X3 X4 xy X2y
1 0.63 1 1 1 0.63 0.63
3 2.05 9 27 81 6.15 18.45
4 4.08 16 64 256 16.32 65.28
6 10.78 36 216 1296 64.68 388.08
∑x=14 ∑y=17.54 ∑x2=62 ∑x3=308 ∑x4=1634 ∑xy=87.78 ∑x2y=472.44
Let,
y = a + bx + cx2
The normal equations are
∑y = ma + b∑x + c∑x 2

∑xy = a∑x + b∑x 2 + 𝑐∑x 3

∑x 2 𝑦 = 𝑎∑x 2 + 𝑏∑x 3 + 𝑐∑x 4

Here m = 4
Then, 4a + 14b + 62c = 17.54
14a + 62b + 308c = 87.78
62a + 308b + 1634c = 472.44
By solving , a = 1.24 , b = -1.05 , c = 0.44

51
The required second degree parabola is,
y = 1.24 − 1.05x + 0.44x 2

Program:

#include<iostream>
#include<iomanip>
#include<cmath>
using namespace std;
int main()
{
cout.precision(6);
cout.setf(ios::fixed);
int i,j,k,n,m;
cout<<"\nEnter the number of points to be entered\n";
cin>>m;
double p[m], q[m], temp, t;
cout<<"\nEnter the x values\n";
for(i=0; i<m; i++)
{
cin>>p[i];
}
cout<<"\nEnter the y values\n";
for(i=0; i<m; i++)
{
cin>>q[i];
}
n=2; //Set degree of polynomial
double a[n+1][n+2],x[n+1];
for(i=0; i<=n; i++)
{
//Calculate different element of the augmented matrix in a row
for(j=0; j<=n; j++)
{
a[i][j]=0;
for(k=0;k<m;k++)
{
a[i][j]=a[i][j]+pow(p[k],(i+j));
}
}
//Calculate last element of the augmented matrix in a row
a[i][n+1]=0;
for(k=0;k<m;k++)
{
a[i][n+1]=a[i][n+1]+pow(p[k],i)*q[k];
}
}
n=n+1;
//Do pivoting
for(k=0; k<n-1; k++)

52
 {
for(i=k+1; i<n; i++)
{
if(fabs(a[k][k])<fabs(a[i][k]))
{
for(j=0; j<=n; j++)
{
temp=a[k][j];
a[k][j]= a[i][j];
a[i][j]=temp;
}
}
}
}
//Do the elementary row operation
for(k=0; k<n-1; k++)
{
for(i=k+1; i<n; i++)
{
t=a[i][k]/a[k][k];
for(j=0; j<=n; j++)
{
a[i][j]=a[i][j]-t*a[k][j];
}
}
}
//Back substitution
for(i=n-1; i>=0; i--)
{
x[i]=a[i][n];
for(j=i+1; j<n; j++)
{
if(j!=i)
x[i]=x[i]-a[i][j]*x[j];
}
x[i]=x[i]/a[i][i];
}
cout<<"\nThe values of the coefficients are\n";
for(i=0; i<n; i++)
{
cout<<x[i]<<endl;
}
cout<<"\nThe required polynomial is\n";
cout<<"y="<<x[0]<<"+("<<x[1]<<")x+"<<x[2]<<"x^2"<<endl;
return 0;
}

Output:
Enter the number of points to be entered
4

53
Enter the x values
1 3 4 6
Enter the y values
0.63 2.05 4.08 10.78
The values of the coefficients are
1.240000
-1.050000
0.440000
The required polynomial is
y= 1.240000 + (-1.050000)x + 0.440000x^2

Problem: 25 7.(c) Curve fitting of the form y = a + bx + cx2+ dx3.

Example: Fit a third degree polynomial to the given data (0, 0), (1, 1.2), (2, 7.3), (3, 25.6),
(4,63.5).
Solution: The table of values is:

x y X2 X3 X4 X5 X6 xy X2y X3y
0 0 0 0 0 0 0 0 0 0
1 1.2 1 1 1 1 1 1.2 1.2 1.2
2 7.3 4 8 16 32 64 14.6 29.2 58.4
3 25.6 9 27 81 243 729 76.8 230.4 691.2
4 63.5 16 64 256 1024 4096 254 1016 4064
10 97.6 30 100 354 1300 4890 346.6 1276.8 4814.8

Let, y = a + bx + cx 2 + dx 3 , here m = 5
And,
∑y = ma + b∑x + c∑x 2 + d∑x 3

∑xy = a∑x + b∑x 2 + 𝑐∑x 3 + 𝑑∑x 4

∑x 2 𝑦 = 𝑎∑x 2 + 𝑏∑x 3 + 𝑐∑x 4 + 𝑑∑x 5

∑x 3 y = a∑x 3 + b∑x 4 + c∑x 5 + d∑x 6

The the normal equations are,

5a + 10b + 30c + 100d = 97.6

10a + 30b + 100c + 354d = 346.6

30a + 100b + 354c + 1300d = 1276.8

54
 100a + 354b + 1300c + 4890d = 4814.8
By solving we get,
𝑎 = −0.001429 , 𝑏 = 1.2178 , 𝑐 = −1.2357 , 𝑑 = 1.2225 ;

Program:

#include<iostream>
#include<iomanip>
#include<cmath>
using namespace std;
int main()
{
cout.precision(6);
cout.setf(ios::fixed);
int i,j,k,n,m;
cout<<"\nEnter the number of points to be entered\n";
cin>>m;
double p[m], q[m], temp, t;
cout<<"\nEnter the x values\n";
for(i=0; i<m; i++)
{
cin>>p[i];
}
cout<<"\nEnter the y values\n";
for(i=0; i<m; i++)
{
cin>>q[i];
}
n=3; //Set degree of polynomial
double a[n+1][n+2],x[n+1];
for(i=0; i<=n; i++)
{
//Calculate different element of the augmented matrix in a row
for(j=0; j<=n; j++)
{
a[i][j]=0;
for(k=0;k<m;k++)
{
a[i][j]=a[i][j]+pow(p[k],(i+j));
}
}
//Calculate last element of the augmented matrix in a row
a[i][n+1]=0;
for(k=0;k<m;k++)
{
a[i][n+1]=a[i][n+1]+pow(p[k],i)*q[k];
}
}
n=n+1;

55
//Do pivoting
for(k=0; k<n-1; k++)
{
for(i=k+1; i<n; i++)
{
if(fabs(a[k][k])<fabs(a[i][k]))
{
for(j=0; j<=n; j++)
{
temp=a[k][j];
a[k][j]= a[i][j];
a[i][j]=temp;
}
}
}
}
//Do the elementary row operation
for(k=0; k<n-1; k++)
{
for(i=k+1; i<n; i++)
{
t=a[i][k]/a[k][k];
for(j=0; j<=n; j++)
{
a[i][j]=a[i][j]-t*a[k][j];
}
}
}
//Back substitution
for(i=n-1; i>=0; i--)
{
x[i]=a[i][n];
for(j=i+1; j<n; j++)
{
if(j!=i)
x[i]=x[i]-a[i][j]*x[j];
}
x[i]=x[i]/a[i][i];
}
cout<<"\nThe values of the coefficients are\n";
for(i=0; i<n; i++)
{
cout<<x[i]<<endl;
}
cout<<"\nThe required polynomial is\n";

cout<<"y="<<x[0]<<"+("<<x[1]<<")x+("<<x[2]<<")x^2+("<<x[3]<<")x^3"<<endl;
return 0;
}

56
Output:
Enter the number of points to be entered
5
Enter the x values
0 1 2 3 4
Enter the y values
0 1.2 7.3 25.6 63.5
The values of the coefficients are
-0.001429
1.217857
-1.235714
1.225000
The required polynomial is
y= -0.001429 + (1.217857)x + (-1.235714)x^2 + (1.225000)x^3
Comments: Since the fitted polynomial obtained by the code is amost same as the fitted polynomial
that we found earlier by manual calculation. So this code could be implement to fitted other
polynomial of degree three for a given set of points.

The End

57