Tutorial 1 : Taylor Series

Problem: Find approximate value of e5 using Taylor series.

Description: The Taylor Series Approximation method is a numerical technique
used to approximate the value of a function at a given point. It is based on the
Taylor series expansion of a function, which represents the function as an infinite
sum of terms.
Method to solve:
Step 1: The General Taylor Series Formula
'' ''' n
f (x 0) f (x 0) f (x 0)
f(x1)=f(x0)+f′(x0)h+h2* +h3* +……..+ hn*
2! 3! n!

Step 2: f(x) = e5
f′(x) = e5
Step 3: Calculate the terms and approximation up to 20 terms.
Step 4: The true value is e5= 148.4131591026.
TV − AV
From that we can find the true relative error by: TV
.

Flow Chart:

Start

Define f(x) and Number

of Terms for calculation.

Find Taylor’s Series,

n k n
f (x−a) TV − AV
e x =∑ ∗(a) and Error =
k=1 n! TV

K>=n
No
Yes
End
Main Program:
clc;
clear all;
x = 5;
num_terms = 20;
taylor_sum = 0;
previous_sum = 0;
convergence = zeros(1, num_terms + 1);
rel_error_plot = zeros(1, num_terms + 1);

for n = 0:num_terms
current_term = (x^n) / gamma(n + 1);
taylor_sum = taylor_sum + current_term;
convergence(n + 1) = taylor_sum;

true_val = exp(x);
rel_error_plot(n + 1) = abs((true_val - taylor_sum) / true_val) * 100;

if n > 0
approx_relative_error = abs((taylor_sum - previous_sum) / taylor_sum);
fprintf('Term %2d: %.10f, Current Sum: %.10f, Abs Rel Error: %.10f (%.2f%
%)\n', ...
n, current_term, taylor_sum, approx_relative_error, approx_relative_error *
100);
else
 fprintf('Term %2d: %.10f, Current Sum: %.10f\n', n, current_term,
taylor_sum);
end
previous_sum = taylor_sum;
end

true_value = exp(x);
true_relative_error = abs((true_value - taylor_sum) / true_value);

fprintf('\n------------------------\nFinal Results:\n------------------------\n');
fprintf('The final approximate value is: %.10f\n', taylor_sum);
fprintf('The true value is: %.10f\n', true_value);
fprintf('The final true relative error is: %.10f (%.2f%%)\n', true_relative_error,
true_relative_error * 100);

figure;
n_axis = 0:num_terms;

[ax, h1, h2] = plotyy(n_axis, convergence, n_axis, rel_error_plot);

hold(ax(1), 'on');
h3 = plot(ax(1), n_axis, true_value * ones(1, num_terms + 1), '--r', 'linewidth', 2);

set(h1, 'LineStyle', '-', 'Marker', 'o', 'LineWidth', 2);

ylabel(ax(1), 'Approximation Value', 'FontSize', 20);
set(h2, 'LineStyle', '-', 'Marker', 'x', 'LineWidth', 1.5, 'Color', [0, 0.5, 0]);
ylabel(ax(2), 'True Relative Error (%)', 'FontSize', 20);
set(ax(2), 'YColor', [0, 0.5, 0]);
xlabel('Number of Terms (n)', 'FontSize', 20);
title(sprintf('Taylor Series Approximation and Error for e^{%.0f}', x), 'FontSize',
20);
grid on;
set(ax(1), 'FontSize', 18);
set(ax(2), 'FontSize', 18);
legend([h1, h3, h2], {'Series Approximation', 'True Value', 'Relative Error
(%)'}, ...
'location', 'northeast', 'FontSize', 18);

hold(ax(1), 'off');
Result:
Term 0: 1.0000000000, Current Sum: 1.0000000000
Term 1: 5.0000000000, Current Sum: 6.0000000000, Abs Rel Error:
0.8333333333 (83.33%)
Term 2: 12.5000000000, Current Sum: 18.5000000000, Abs Rel Error:
0.6756756757 (67.57%)
Term 3: 20.8333333333, Current Sum: 39.3333333333, Abs Rel Error:
0.5296610169 (52.97%)
Term 4: 26.0416666667, Current Sum: 65.3750000000, Abs Rel Error:
0.3983428936 (39.83%)
Term 5: 26.0416666667, Current Sum: 91.4166666667, Abs Rel Error:
0.2848678213 (28.49%)
Term 6: 21.7013888889, Current Sum: 113.1180555556, Abs Rel Error:
0.1918472589 (19.18%)
Term 7: 15.5009920635, Current Sum: 128.6190476190, Abs Rel Error:
0.1205186351 (12.05%)
Term 8: 9.6881200397, Current Sum: 138.3071676587, Abs Rel Error:
0.0700478522 (7.00%)
Term 9: 5.3822889109, Current Sum: 143.6894565697, Abs Rel Error:
0.0374577860 (3.75%)
Term 10: 2.6911444555, Current Sum: 146.3806010251, Abs Rel Error:
0.0183845703 (1.84%)
Term 11: 1.2232474798, Current Sum: 147.6038485049, Abs Rel Error:
0.0082873685 (0.83%)
Term 12: 0.5096864499, Current Sum: 148.1135349548, Abs Rel Error:
0.0034411875 (0.34%)
Term 13: 0.1960332500, Current Sum: 148.3095682048, Abs Rel Error:
0.0013217842 (0.13%)
Term 14: 0.0700118750, Current Sum: 148.3795800797, Abs Rel Error:
0.0004718431 (0.05%)
Term 15: 0.0233372917, Current Sum: 148.4029173714, Abs Rel Error:
0.0001572563 (0.02%)
Term 16: 0.0072929036, Current Sum: 148.4102102750, Abs Rel Error:
0.0000491402 (0.00%)
Term 17: 0.0021449717, Current Sum: 148.4123552467, Abs Rel Error:
0.0000144528 (0.00%)
Term 18: 0.0005958255, Current Sum: 148.4129510722, Abs Rel Error:
0.0000040146 (0.00%)
Term 19: 0.0001567962, Current Sum: 148.4131078683, Abs Rel Error:
0.0000010565 (0.00%)
Term 20: 0.0000391990, Current Sum: 148.4131470674, Abs Rel Error:
0.0000002641 (0.00%)
------------------------
Final Results:
------------------------
The final approximate value is: 148.4131470674
The true value is: 148.4131591026
The final true relative error is: 0.0000000811 (0.00%)

*Chart is attached separately.

Conclusion: The Taylor series approximation of e 5 demonstrates that numerical
accuracy increases significantly as more terms are added to the expansion,
transitioning from an initial relative error of over 80% to nearly 0.00001% within
20 iterations.
 Tutorial 2 : Bisection Method
Problem-I: Find approximate root using Bisection method for f(x)=x5-5x-1 up to
error 10-6.
Description: The Bisection method is a numerical technique for finding the roots
of a continuous function. It works by repeatedly halving the interval between two
points, one with a positive function value and the other with a negative one.
Method to solve:
Step 1: Define f(x), first approximate value, second approximate value, True value
and error.
Step 2: Check the value of f(xi)*f(xf). If value is greater than zero then display
error and we have to change our approximated value. If the value is less than zero
then find the value of mean of x i and xf. Also find the error. If error is less than the
assigned value then we have to stop. Store the value of x m, error and number of
iteration to display in output and graph.
Step 3: Print the value of xm and error.
Step 4: Plot the graph of xm v/s iteration and error v/s iteration.
Flowchart: - START

Define f(x)
Initialize a and b
Solution does
not lie in
between a &
[Link]
No value of a & b
f(a)f(b)<0

yes

C = (a + b)/2

no
f(a)f(c)<0
c=a

yes

c=b

no
error<Tol

yes

STOP
Main Program:
%% Bisection Method for x^5 - 5x - 1 = 0
clear all
clc
f = @(x) x.^5 - 5.*x - 1;
a = 1;
b = 2;
tol = 1e-6;
max_iter = 100;
iter_vals = [];
c_vals = [];
relative_errors = [];
c_old = 0;
if f(a) * f(b) > 0
fprintf('The initial interval [a, b] does not contain a root.\n');
else
fprintf('Iter. No. Approximate Value Absolute Relative Error\n');
for i = 1:max_iter
c = (a + b) / 2;
if i == 1
rel_error = NaN;
else
rel_error = abs((c - c_old) / c);
end
c_old = c;
 fprintf('%-10d %-25.6f %-25.6f\n', i, c, rel_error);
iter_vals = [iter_vals, i];
c_vals = [c_vals, c];
relative_errors = [relative_errors, rel_error];
if rel_error < tol && i > 1
fprintf('Root found after %d iterations.\n', i);
fprintf('Final Approximate Root: %.6f\n', c);
break;
end

% Update the interval for the next iteration

if f(c) * f(a) < 0
b = c;
else
a = c;
end
end
end

figure;
subplot(2,1,1);
plot(iter_vals, c_vals, 'b-o', 'LineWidth', 1.5, 'MarkerSize', 6);
grid on;
title('Approximate Value Convergence', 'fontsize', 20);
xlabel('Number of Iterations', 'fontsize', 20);
ylabel('Approximate Value', 'fontsize', 20);
subplot(2,1,2);
semilogy(iter_vals(2:end), relative_errors(2:end), 'r-o', 'LineWidth', 1.5,
'MarkerSize', 6);
grid on;
title('Absolute Relative Error Convergence', 'fontsize', 20);
xlabel('Number of Iterations', 'fontsize', 20);
ylabel('Absolute Relative Error', 'fontsize', 20);

Result:
Iter. No. Approximate Value Absolute Relative Error
1 1.500000 NaN
2 1.750000 0.142857
3 1.625000 0.076923
4 1.562500 0.040000
5 1.531250 0.020408
6 1.546875 0.010101
7 1.539063 0.005076
8 1.542969 0.002532
9 1.541016 0.001267
10 1.541992 0.000633
11 1.541504 0.000317
12 1.541748 0.000158
13 1.541626 0.000079
14 1.541687 0.000040
15 1.541656 0.000020
16 1.541641 0.000010
17 1.541649 0.000005
18 1.541653 0.000002
19 1.541651 0.000001
20 1.541652 0.000001
Root found after 20 iterations.
Final Approximate Root: 1.541652

*Chart is attached separately.

Conclusion: The Bisection Method successfully identified the root of the equation
x5 - 5x - 1 = 0 within the interval [1, 2], demonstrating consistent convergence
toward the final approximate value..
Problem-II: Find approximate root using Regula-Falsi / False position / Linear
interpolation method for f(x)=x5-5 up to error 10-6.
Description: The Regula Falsi method starts with an initial interval [x i,xf] where
the function values, f(xi) and f(xf), have opposite signs. This ensures that a root
exists within the interval. Unlike the bisection method, which simply halves the
interval, the Regula Falsi method uses a line connecting the two points (x i,f(xi))
and (xf,f(xf)). The point where this line intersects the x-axis is calculated as the new
approximation of the root, xr.
Method to solve:
Step 1: Define f(x), first approximate value, second approximate value, True value
and error.
Step 2: Check the value of f(xi)*f(xf). If value is greater than zero then display
error and we have to change our approximated value. If the value is less than zero
then find the value of next approximate value by f(x i) and f(xf). Also find the error.
If error is less than the assigned value then we have to stop. Store the value of x r,
error and number of iteration to display in output and graph.
Step 3: Print the value of xr and error.
Step 4: Plot the graph of xr v/s iteration and error v/s iteration.
Flowchart: - START

Define f(x)
Solution does
Initialize x0 and x1 not lie in
between x0 &
x1.
Change value
No of x0 & x1
f(x0)f(x1)<0

Yes

x 0 f (x 1) – x 1 f ( x 0)
x 2=
x 1+ x 0

No
f(x2)f(x1)<0
x1 = x0 & x0
= x2

Yes

x0 = x1 & x1 = x2

No
error<Tol

Yes

STOP
Main Program:
clc;
clear all;
f = @(x) x.^5 - 5;
xi = 1.0;
xf = 2.0;
xs = 5^(1/5);
e = 10^-7;
n = 100;
if (f(xi) * f(xf) < 0)
vxr = zeros(1, n);
verr = zeros(1, n);

for i = 1:n
xr = xf - (f(xf) * (xi - xf)) / (f(xi) - f(xf));
vxr(i) = xr;

errt = 100 * abs((xs - xr) / xs);

verr(i) = errt;

fprintf('ITR: %d, x = %0.8f, True Relative Error = %0.8f%%\n', i, xr, errt);

if errt < e
break;
end
 if f(xi) * f(xr) < 0
xf = xr;
else
xi = xr;
end
end

fprintf('\nFinal root found at x = %0.8f\n', xr);

fprintf('Number of iterations: %d\n', i);

figure;
subplot(2, 1, 1);
plot(1:i, vxr(1:i), 'b-o', 'LineWidth', 1.5);
xlabel('Number of Iterations', 'FontSize', 18);
ylabel('Value (x)', 'FontSize', 18);
title('Regula Falsi Method Convergence', 'FontSize', 20);
grid on;
subplot(2, 1, 2);
semilogy(1:i, verr(1:i), 'r-x', 'LineWidth', 1.5);
xlabel('Number of Iterations', 'FontSize', 18);
ylabel('Error (%)', 'FontSize', 18);
title('True Relative Error (%)', 'FontSize', 20);
grid on;
else
fprintf('Error: f(xi) and f(xf) must have opposite signs.\n');
end
Result:
ITR: 1, x = 1.12903226, True Relative Error = 18.17003797%
ITR: 2, x = 1.22042813, True Relative Error = 11.54585106%
ITR: 3, x = 1.28144039, True Relative Error = 7.12380623%
ITR: 4, x = 1.32032415, True Relative Error = 4.30559101%
ITR: 5, x = 1.34430821, True Relative Error = 2.56727500%
ITR: 6, x = 1.35878772, True Relative Error = 1.51782951%
ITR: 7, x = 1.36741227, True Relative Error = 0.89273914%
ITR: 8, x = 1.37250741, True Relative Error = 0.52345446%
ITR: 9, x = 1.37550270, True Relative Error = 0.30636183%
ITR: 10, x = 1.37725843, True Relative Error = 0.17911009%
ITR: 11, x = 1.37828581, True Relative Error = 0.10464764%
ITR: 12, x = 1.37888638, True Relative Error = 0.06111913%
ITR: 13, x = 1.37923725, True Relative Error = 0.03568867%
ITR: 14, x = 1.37944217, True Relative Error = 0.02083667%
ITR: 15, x = 1.37956182, True Relative Error = 0.01216450%
ITR: 16, x = 1.37963168, True Relative Error = 0.00710135%
ITR: 17, x = 1.37967246, True Relative Error = 0.00414550%
ITR: 18, x = 1.37969627, True Relative Error = 0.00241995%
ITR: 19, x = 1.37971017, True Relative Error = 0.00141264%
ITR: 20, x = 1.37971828, True Relative Error = 0.00082462%
ITR: 21, x = 1.37972302, True Relative Error = 0.00048137%
ITR: 22, x = 1.37972578, True Relative Error = 0.00028100%
ITR: 23, x = 1.37972740, True Relative Error = 0.00016403%
ITR: 24, x = 1.37972834, True Relative Error = 0.00009575%
ITR: 25, x = 1.37972889, True Relative Error = 0.00005589%
ITR: 26, x = 1.37972921, True Relative Error = 0.00003263%
ITR: 27, x = 1.37972940, True Relative Error = 0.00001905%
ITR: 28, x = 1.37972951, True Relative Error = 0.00001112%
ITR: 29, x = 1.37972957, True Relative Error = 0.00000649%
ITR: 30, x = 1.37972961, True Relative Error = 0.00000379%
ITR: 31, x = 1.37972963, True Relative Error = 0.00000221%
ITR: 32, x = 1.37972964, True Relative Error = 0.00000129%
ITR: 33, x = 1.37972965, True Relative Error = 0.00000075%
ITR: 34, x = 1.37972966, True Relative Error = 0.00000044%
ITR: 35, x = 1.37972966, True Relative Error = 0.00000026%
ITR: 36, x = 1.37972966, True Relative Error = 0.00000015%
ITR: 37, x = 1.37972966, True Relative Error = 0.00000009%

Final root found at x = 1.37972966

Number of iterations: 37

*Chart is attached separately.

Conclusion: The Regula Falsi method demonstrates an efficient approach to root-
finding by utilizing linear interpolation between interval boundaries, often
resulting in faster convergence than the standard Bisection method.
By incorporating the function's magnitude at each endpoint, the algorithm
intelligently narrows the search space toward the actual root of x 5-5=0. The
simulation results confirm that the true relative error decreases consistently,
reaching the desired precision of 10-7 within a finite number of iterations.
 Tutorial 3: Gauss Elimination Method
Problem: solve the following system of equations by gauss elimination method:
2x + y - z = 8
-3x - y + 2z = -11
-2x + y + 2z = -3
Flow Chart:

Start

Define Matrix

Make lower triangular matrix

zero by row operations

Do backward matrix
multiplication process and
find all unknowns

Print Results

End
Description: The system of equations is first represented as an augmented matrix,
which combines the coefficients of the variables and the constants. Using a series
of elementary row operations an upper triangular matrix with zeros below the main
diagonal.
Method to solve:
Step 1: Form the augmented matrix.
Step 2: Take elementary row operation and make lower triangle zero.
Step 3: Solve the equation.
Main Program:
%% Gauss Elimination Method for a 3-variable system
clc;
clear all;
A = [2, 1, -1;
-3, -1, 2;
-2, 1, 2];
b = [8; -11; -3];
M = [A, b];
n = size(A, 1);
for k = 1:n-1
[~, p] = max(abs(M(k:n, k)));
p = p + k - 1;

if p ~= k
M([k, p], :) = M([p, k], :);
end
 for i = k+1:n
m = M(i, k) / M(k, k);
M(i, k:n+1) = M(i, k:n+1) - m * M(k, k:n+1);
end
end

x = zeros(n, 1);
x(n) = M(n, n+1) / M(n, n);

for i = n-1:-1:1
x(i) = (M(i, n+1) - M(i, i+1:n) * x(i+1:n)) / M(i, i);
end

fprintf('The augmented matrix after forward elimination is:\n');

disp(M);

fprintf('The solution to the system of equations is:\n');

labels = ['x', 'y', 'z'];
for i = 1:n
fprintf('%c = %.4f\n', labels(i), x(i));
end

Result:
The augmented matrix after forward elimination is:
-3.0000 -1.0000 2.0000 -11.0000
 0 1.6667 0.6667 4.3333
0 0 0.2000 -0.2000

The solution to the system of equations is:

x = 2.0000
y = 3.0000
z = -1.0000

Conclusion:
The Gaussian elimination method with partial pivoting provides a robust numerical
approach for solving linear systems by transforming an augmented matrix into an
upper triangular form. By incorporating row-swapping logic to select the largest
absolute pivot element, the algorithm effectively prevents division by zero and
significantly minimizes numerical round-off errors.
 Tutorial 4 : Newton Raphson Method
Problem: solve the approximate root of equation by Newton Raphson method:
f(x)=x5−5x−5.

Flow Chart:

Start

Define function and its

differentiation, initialize
error, initial value

Calculate next
approximate value

No Error<Es

Yes

Print Root and error

End
Description: The Newton-Raphson method is a powerful and widely used iterative
technique for finding the roots (or zeros) of a non-linear equation, f(x)=0. The
method starts with an initial guess and uses the tangent line to the function at that
point to find a better, closer approximation of the root.

Method to solve:

The Newton-Raphson method is an efficient iterative technique used to find the

roots of a real-valued function f(x). The iterative formula is:

Xn+1=xn− f(xn)/ f′(xn)

The process implemented in the code is as follows:

1. Define Function and Derivative: The function f(x)=x5−5x−5 and its

derivative f′(x)=5x4−5 are defined.
2. Initialize Parameters: An initial guess of a=10 and an error tolerance (tol) of
10−5 are set.
3. Iteration Loop: The process enters a loop that continues as long as the
approximate relative error (ea) is greater than the tolerance.
4. Calculate New Root: The new approximation (b) is calculated using the
Newton-Raphson formula.
5. Calculate Error: The approximate relative error is calculated as e a=100×
|b−ab|.
6. Store and Update: The iteration number, error, and new root are stored, and
the old root (a) is updated to the new root (b) for the next iteration.
7. Print Results: The final approximate root and the required number of
iterations are printed upon convergence.

Main Program:
% Newton Raphson method
clc; clear all;
f= @(x) x^5-5*x-5;
fp =@(x) 5*x^4-5;
a=10;
i=0;
ea=100;
tol=1e-5;
iter=[];
error=[];
root=[];

while ea > tol

i = i + 1;
b = a - (f(a) / fp(a));
if b ~= 0
ea = 100 * abs((b - a) / b);
else
ea = 0;
end
iter(i)=i;
error(i)=ea;
root(i)=b;
a = b;
fprintf('At iteration %d\t approximate root is %12.8f\t with error %12.6f\n', i, b,
ea);
end
fprintf('\n Approximate root of the given equation is %f at iteration %d \n', b,i);

figure;
% Subplot 1: Iteration vs Root
subplot(2,1,1);
plot(iter,root,'-o','color','r', 'LineWidth', 1.5);
title('Newton Raphson Method:Iteration vs Root', 'FontSize', 20);
xlabel('Iteration', 'FontSize', 18);
ylabel('Root Value', 'FontSize', 18);
set(gca, 'FontSize', 14);
grid on;
% Subplot 2: Iteration vs Error
subplot(2,1,2);
plot(iter,error,'-o', 'LineWidth', 1.5);
title('Newton Raphson Method:Iteration vs Error', 'FontSize', 20);
xlabel('Iteration', 'FontSize', 18);
ylabel('Error (%)', 'FontSize', 18);
set(gca, 'FontSize', 14);
grid on;

Result:
At iteration 1 approximate root is 8.00090009 with error 24.985938
At iteration 2 approximate root is 6.40252652 with error 24.964732
At iteration 3 approximate root is 5.12566664 with error 24.911099
At iteration 4 approximate root is 4.10793351 with error 24.774820
At iteration 5 approximate root is 3.30145184 with error 24.428092
At iteration 6 approximate root is 2.67207091 with error 23.554050
At iteration 7 approximate root is 2.20043614 with error 21.433695
At iteration 8 approximate root is 1.88333620 with error 16.837140
At iteration 9 approximate root is 1.72311766 with error 9.298177
At iteration 10 approximate root is 1.68281441 with error 2.394991
At iteration 11 approximate root is 1.68050131 with error 0.137643
At iteration 12 approximate root is 1.68049401 with error 0.000434
At iteration 13 approximate root is 1.68049401 with error 0.000000

Approximate root of the given equation is 1.680494 at iteration 13

*Chart is attached separately.

Conclusion: The Newton-Raphson method successfully approximates

the root of the nonlinear equation x5 - 5x - 5 = 0 by utilizing the
function's derivative to achieve rapid convergence. The simulation
demonstrates the method's characteristic quadratic convergence, as the
approximate relative percentage error drops significantly within a few
iterations once the value approaches the actual root.
 Tutorial 5 : Secant Method
Problem: solve the approximate root of equation by secant method: f(x) =e− x −x.
Flow Chart:

Start

Define function,
initialize error, initial
value

Calculate next
approximate value

No
Error<εs

Yes

Print Root and error

End
Description: When the derivatives of the function is not known, we can use secant method
instead of Newton Raphson method. We have to just put derivative part in NR method by
forward, backward or central difference. Here, we have used backward difference.

Method to solve:
The Secant Method can be viewed as an approximation of the Newton-Raphson
method where the derivative is replaced by a finite-difference approximation. The
iterative formula is:
f (b)∗(b−a)
c= b - ( f (b)−f (a))

The process implemented in the code is as follows:

1. Define Function: The function f(x)=e− x −x.
2. Initialize Parameters: An initial guess of a=1, b=2 and an error tolerance
(tol) of 10−5 are set.
3. Iteration Loop: The process enters a loop that continues as long as the
approximate relative error (ea) is greater than the tolerance.
4. Calculate New Root: The new approximation (c) is calculated using the
secant formula.
5. Calculate Error: The approximate relative error is calculated as

ea=100× b|b−a|
6. Store and Update: The iteration number, error, and new root are stored, and
the old root (a) is updated to the new root (b) for the next iteration.
7. Print Results: The final approximate root and the required number of
iterations are printed upon convergence
Main Program:
% Tutorial 5
% Secant method
clc; clear all;
f= @(x) exp(-x)-x;
a=1;
b=2;
i=0;
ea=100;
tol=1e-5;
iter=[];
error=[];
root=[];

while ea > tol

i = i + 1;
c = b - (f(b)*(b-a)) / (f(b)-f(a));
if c ~= 0
ea = 100 * abs(c - b) / abs(c);
end
a=b;
b=c;
iter(i)=i;
error(i)=ea;
root(i)=c;
fprintf('At iteration %d\t approximate root is %12.8f\t with error %12.6f\n', i, c,
ea);
end

fprintf('\n Approximate root of the given equation is %f at iteration %d \n', c,i);

figure;
% Subplot 1: Root Convergence
subplot(2,1,1);
plot(iter,root,'-o','color','r', 'LineWidth', 1.5);
title('Secant Method: Iteration vs Root', 'FontSize', 20);
xlabel('Iteration Number', 'FontSize', 18);
ylabel('Root Value', 'FontSize', 18);
grid on;

% Subplot 2: Error Convergence

subplot(2,1,2);
plot(iter,error,'-o', 'LineWidth', 1.5);
title('Secant Method: Iteration vs Error', 'FontSize', 20);
xlabel('Iteration Number', 'FontSize', 18);
ylabel('Error (%)', 'FontSize', 18);
grid on;

Result:
At iteration 1 approximate root is 0.48714165 with error 310.558199
At iteration 2 approximate root is 0.58377969 with error 16.553853
At iteration 3 approximate root is 0.56738645 with error 2.889254
At iteration 4 approximate root is 0.56714256 with error 0.043003
At iteration 5 approximate root is 0.56714329 with error 0.000129
At iteration 6 approximate root is 0.56714329 with error 0.000000
Approximate root of the given equation is 0.567143 at iteration 6.
Conclusion:
The Secant method effectively approximates the root of the transcendental
equation e-x - x = 0 by using a finite difference approximation of the derivative
through two initial guesses. The simulation results demonstrate a super-linear
convergence rate, where the approximate relative error is reduced to less than 10 -5
within a few iterations without requiring the analytical derivative of the function.
 Tutorial 6 : Newton’s Divided Difference Interpolation Method

Problem: Estimate the value of log(2.5) using Newton’s Divided Difference

Interpolation Method, and determine the true error.

Flow Chart:

Start

Define function,
initialize variables and
True value.

Calculate divided differences and store

value in b0,b1,b2…bn

Calculate approximate value and true

error at each iteration

Print solution and error

End
Description: Interpolation is the process of estimating the value of a function for any
intermediate value of the independent variable when the function is known at a set of discrete
points.

Newton’s Divided Difference Method is a form of polynomial interpolation that can be used for
unequally spaced data points.

Method to solve:
Input data points (x0,f(x0)),(x1,f(x1)),…,(xn,f(xn )) and the value xn where
interpolation is required.
Initialize a divided difference table dd(n,n) and store the first column as the given
f(xi)f(x_i)f(xi) values.
Compute higher-order divided differences using the recursive formula.
Construct the interpolation polynomial incrementally using:
P=dd(1,1)+dd(1,2)(x−x0)+dd(1,3)(x−x0)(x−x1)+...P = dd(1,1) + dd(1,2)(x - x_0) +
dd(1,3)(x - x_0)(x - x_1) + ...P=dd(1,1)+dd(1,2)(x−x0)+dd(1,3)(x−x0)(x−x1)+...
Evaluate P(xn)P(x_n)P(xn) for the required point.
Compare with the true value (if known) to compute the error.
Main Program:
% Newton Divided Difference Interpolation
clc;
clear;
close all;
x = [1 2 3 4 5 6 7];
y = log(x);
xn = 2.5;
n = length(x);
dd = zeros(n, n);
dd(:, 1) = y';
for j = 2:n
for i = 1:n - j + 1
dd(i, j) = (dd(i+1, j-1) - dd(i, j-1)) / (x(i+j-1) - x(i));
end
end

TrueValue = log(xn);
P = dd(1, 1);
product_term = 1;
ApproximateValues = zeros(1, n);
TrueError = zeros(1, n);
ApproximateValues(1) = P;
TrueError(1) = abs(TrueValue - P);
for j = 2:n
product_term = product_term * (xn - x(j - 1));
P = P + dd(1, j) * product_term;
ApproximateValues(j) = P;
TrueError(j) = abs(TrueValue - P);
end

for j = 1:n
k = j - 1;
fprintf('For order %d Approximate solution is %0.10f with true error %0.10e\n',
k, ApproximateValues(j), TrueError(j));
end
disp(['True value log(', num2str(xn), '): ', num2str(TrueValue)]);
disp(['Final Interpolated Value P(', num2str(xn), '): ', num2str(P)]);
disp(['Final True Error: ', num2str(TrueError(n))]);
disp(' ');

Result:
For order 0 approximate solution is 0.000000000 with true error 9.1629073187e-
01
For order 1 approximate solution is 1.039720770 with true error 1.2343003897e-
01
For order 2 approximate solution is 0.931839993 with true error 1.5549261796e-
02
For order 3 approximate solution is 0.921221303 with true error 4.9305719966e-
03
For order 4 approximate solution is 0.918487213 with true error 2.1964816308e-
03
For order 5 approximate solution is 0.917466199 with true error 1.1754675561e-
03
For order 6 approximate solution is 0.916996816 with true error 7.0608504751e-
04
True value log(2.5): 0.91629
Final Interpolated Value P(2.5): 0.917
Final True Error: 0.00070609

Conclusion:
The Newton’s Divided Difference Method effectively estimates intermediate
function values even for unequally spaced data. The accuracy improves with more
data points and depends on the smoothness of the function.
 Tutorial 7 : Lagrange Interpolation

Problem: Estimate the value of log (2.5) using the Lagrange Interpolation Method,
and determine the final true relative percentage error.
The data points used are:
x = [1, 2, 3, 4, 5, 6, 7]
f(x)=log(x)

Flow Chart:

Start

Define function, initialize

variables and True value.

Calculate Lagrange basis polynomial

n
x−xj
Li(x)= ∏ xi−xj
j =1 , j ≠1

Calculate approximate value and true

error at each iteration
n
Pn-1(x )=∑ yi∗Li (x )
i=1

Print solution and error

End
Description: The Lagrange Interpolation Method constructs a single polynomial,
Pn−1(x), of degree n−1 that passes exactly through all n given data points. This
polynomial is unique and can be used to approximate the function's value at any
point xn.
Method to solve:
The Lagrange Interpolating Polynomial, Pn−1(x), is defined by the formula:
n

Pn-1(x)=∑ nyi∗Li (x)

i=1

Where Li(x) is the Lagrange basis polynomial, given by:

n
x−xj
Li(x)= ∏ xi−xj
j =1 , j ≠1

 Input Data: Define the input nodes x i and the corresponding function values y i
=log(xi). Define the interpolation point xn=2.5.
 Calculate Basis Polynomials: For each data point x i, calculate the scalar value of
the Lagrange basis polynomial Li(xn).
 Sum the Terms: Calculate the product yi⋅Li(xn) for each i, and sum these terms
to find the final interpolated value, AV=Pn−1(xn).
 Compute Error: Calculate the true value, TV=log(x n), and compute the True
Relative Percentage Error:
TV − AV
Error=¿ TV
∨¿*100

Main Program:
clc;
clear;
close all;
x = [1 2 3 4 5 6 7];
y = log(x);
xn = 2.5;
Li = 1;
n = length(x);
AV = 0;
tv = log(xn);
for i = 1:n
Li = 1;
for j = 1:n
if i ~= j
L_fraction = (xn - x(j)) / (x(i) - x(j));
Li = Li * L_fraction;
end
end
term = y(i) * Li;
AV = AV + term;
error = abs((tv - AV) / tv) * 100;
fprintf('For order %d Approximate solution is %.8f with error %.8f\n',i,AV,tv)
end

fprintf('\n--- Lagrange Interpolation Results ---\n');

fprintf('Interpolation Point (xn): %.2f\n', xn);
fprintf('True Value (log(%.2f)): %.8f\n', xn, tv);
fprintf('Approximate Value (AV): %.8f\n', AV);
fprintf('True Relative Percentage Error: %.5f%%\n', error);
Result:
For order 1 Approximate solution is 0.00000000 with error 0.91629073
For order 2 Approximate solution is 0.25586878 with error 0.91629073
For order 3 Approximate solution is 1.26972485 with error 0.91629073
For order 4 Approximate solution is 0.70112755 with error 0.91629073
For order 5 Approximate solution is 0.99818201 with error 0.91629073
For order 6 Approximate solution is 0.90369470 with error 0.91629073
For order 7 Approximate solution is 0.91699682 with error 0.91629073
--- Lagrange Interpolation Results ---
Interpolation Point (xn): 2.50
True Value (log(2.50)): 0.91629073
Approximate Value (AV): 0.91699682
True Relative Percentage Error: 0.07706%
Conclusion:
The Lagrange Interpolation Method successfully estimated the value of log (2.5)
using the given data points.
The calculated value of log (2.5) using the 6th-order polynomial (since n=7 points)
is 0.91699682.
The True Relative Percentage Error is very small, measuring 0.07706%, indicating
a highly accurate approximation of the log function at the point x=2.5.
 Tutorial 8: Refrigerant Properties

Problem: The objective is to determine how the Coefficient of Performance (COP)

of a standard vapor compression refrigeration cycle changes when the condensing
temperature (T3) is varied, while keeping the evaporating pressure (P1),
condensing pressure (P2), and the degree of superheat (DOS) constant. This
comparative analysis is performed for four common refrigerants: R134A,
Ammonia (R717), Propane (R290), and R123.
Flow Chart:

Start

Initialize different fluids, both

pressure, and DOS

Get different data like T,S,H from REFPROP

Initialize the range for the atmospheric

temperature (T3)

Calculate COP and plot

graph

End
Method to solve:

 Define different fluids and initial pressure and temperature.

 Get the value of T1_Sat at given pressure P1 from refprop.
 Find temperature T1=T1_sat+Degree of Superheat.
 Find the value of Enthalpy and entropy at pressure P1 and temperature T1.
 Process 1-2 is Isentropic, So S1=S2.
 Find the H2 and T2 from pressure P2 and S2.
 3-4 is isenthalpic process, So H3=H4, Find H3 from P2 and T3.
h 1−h 4
 Calculate COP = h 2−h1
 Plot the graph.

Main Program:

clc;
clear all;
Refrigerants = { 'R134A', 'AMMONIA', 'PROPANE', 'R123'};
Fluid_Names = {'R134A', 'R717 (Ammonia)', 'R290 (Propane)', 'R123'};

P1_bar = 1.2;
DOS = 20;
P2_bar = 15;

P1 = P1_bar * 100;
P2 = P2_bar * 100;
P3 = P2;
P4 = P1;

T3_range = 278:1:323;
T3_C_range = T3_range - 273.15;

figure;
hold on;
Colors = lines(length(Refrigerants));

for fluid_index = 1:length(Refrigerants)

Fluid = Refrigerants{fluid_index};
Legend_Name = Fluid_Names{fluid_index};
COP_values = zeros(size(T3_range));

disp(['--- Calculating for ', Legend_Name, ' (', Fluid, ') ---']);

T_sat = refpropm('T', 'P', P1, 'Q', 0, Fluid);

T1 = T_sat + DOS;

H1 = refpropm('H', 'T', T1, 'P', P1, Fluid);

S1 = refpropm('S', 'T', T1, 'P', P1, Fluid);

S2 = S1;
H2 = refpropm('H', 'P', P2, 'S', S2, Fluid);
Q_comp = H2 - H1;

for i = 1:length(T3_range)
T3 = T3_range(i);
 H3 = refpropm('H', 'T', T3, 'P', P3, Fluid);
H4 = H3;
COP = (H1 - H4) / Q_comp;
COP_values(i) = COP;
end

plot(T3_C_range, COP_values, 'Color', Colors(fluid_index,:), 'LineWidth', 2,

'DisplayName', Legend_Name);
end
hold off;
xlabel('Condensing Temperature T3 (°C)');
ylabel('Coefficient of Performance (COP)');
title(['COP vs. Condensing Temperature for Selected Refrigerants']);
legend show;
grid on;
*Chart is attached separately.
Conclusion:

R123 exhibits the highest COP across the entire operating range, starting near 4.8
and consistently remaining superior to the other fluids. This indicates that for this
specific set of operating pressures (1.2 bar evaporating, 15 bar condensing), R123
is the most energy-efficient choice.
R134A and R717 (Ammonia) demonstrate competitive performance in the low to
mid-range temperatures (COP between 2.8 and 3.4). R290 (Propane) shows the
lowest practical COP among the group at lower temperatures, hovering around 2.7.
The selection of the optimal refrigerant depends heavily on the expected
condensing temperature. While R123 offers the maximum energy efficiency under
moderate conditions, its performance stability is limited to condensing
temperatures below approximately 48 °C.
 Tutorial 9: First order and Second Order Regression

Problem: Find the best fit line curve and coefficient of correlation for following
data for first order and second order. Also compare charts.
X1 0 2 2.5 1 4 7
X2 0 1 2 3 6 2
Yi 5 10 9 0 3 27

Description: Curve Fitting is a numerical technique used to construct a

mathematical model that best fits a set of data points. The goal is to find a
continuous function that approximates the data points, allowing for interpolation,
extrapolation, and prediction.

Flow Chart:

START

Input X1 X2 Yi
Setup normal equations
for quadratic fit

Calculate coefficients a0 a1 a2
Y = a0 + a1x + a2x2

Result – find correlation coefficient

(r2)
 Input X1 X2 Yi
Setup normal equations
for linear fit

Calculate coefficients a0 a1
Y = a0 + a1x

Result – find correlation coefficient

(r1)

Compare r2 and r1 and plot graph

STOP
Main Program:
clc;
clear all;
x1 = [0; 2; 2.5; 1; 4; 7];
x2 = [0; 1; 2; 3; 6; 2];
y = [5; 10; 9; 0; 3; 27];

A = [6,sum(x1), sum(x2), sum(x1.^2), sum(x2.^2),sum(y);

sum(x1), sum(x1.^2), sum(x1.*x2), sum(x1.^3),sum(x2.^2.*x1), sum(x1.*y);
sum(x2),sum(x1.*x2),sum(x2.^2),sum((x1.^2).*x2),sum(x2.^3), sum(x2.*y) ;
sum(x1.^2),sum(x1.^3),sum((x1.^2).*x2),sum(x1.^4),sum((x1.^2).*(x2.^2)),sum
((x1.^2).*y);
sum(x2.^2),sum((x1.^2).*x2),sum(x2.^3),sum((x1.^2).*(x2.^2)),sum(x2.^4),sum
((x2.^2).*y)];

n = size(A, 1);

for k = 1:n-1
[~, max_row] = max(abs(A(k:n, k)));
max_row = max_row + k - 1;

A([k, max_row], :) = A([max_row, k], :);

if A(k, k) == 0
error('Matrix is singular or ill-conditioned. Cannot proceed.');
end
 for i = k+1:n
factor = A(i, k) / A(k, k);
A(i, k:n+1) = A(i, k:n+1) - factor * A(k, k:n+1);
end
end

a = zeros(n, 1);

for i = n:-1:1
sum_term = 0;
for j = i+1:n
sum_term = sum_term + A(i, j) * a(j);
end

a(i) = (A(i, n+1) - sum_term) / A(i, i);

end

fprintf('A Second order fit line equation is :\n Yi=%f + %f X1+ %f X2+%f X1^2 +
%f X2^2\n',a(1),a(2),a(3),a(4),a(5))
Yi=a(1)+a(2).*x1+a(3).*x2+a(4).*x1.^2+a(5).*x2.^2;

B= [1,sum(x1),sum(x2),sum(y);
sum(x1), sum(x1.^2), sum(x1.*x2),sum(x1.*y);
sum(x2),sum(x1.*x2), sum(x2.^2),sum(x2.*y)];
n = size(B, 1);
for k = 1:n-1
[~, max_row] = max(abs(B(k:n, k)));
max_row = max_row + k - 1;

B([k, max_row], :) = B([max_row, k], :);

if B(k, k) == 0
error('Matrix is singular or ill-conditioned. Cannot proceed.');
end

for i = k+1:n
factor = B(i, k) / B(k, k);
B(i, k:n+1) = B(i, k:n+1) - factor * B(k, k:n+1);
end
end

b = zeros(n, 1);

for i = n:-1:1
sum_term = 0;
for j = i+1:n
sum_term = sum_term + B(i, j) * b(j);
end
 b(i) = (B(i, n+1) - sum_term) / B(i, i);
end

fprintf('A first order best fit line equation is :\n Zi=%f + %f X1+ %f X2\
n',b(1),b(2),b(3))
Zi=b(1)+b(2).*x1+b(3).*x2;

n=size(A,1);

y_bar=(sum(y))/n;
St1=sum((Zi-y_bar).^2);
Sr1=sum((y-Zi).^2);
r1=((St1-Sr1)/St1)^0.5;
fprintf('\n Coefficient of Corelation is :%f\n', r1);
St2=sum((Yi-y_bar).^2);
Sr2=sum((y-Yi).^2);
r2=((St2-Sr2)/St2)^0.5;
fprintf('\n Coefficient of Corelation is :%f\n', r2);
n=[1,2,3,4,5,6];
plot(n,Yi,'m-p', 'DisplayName', 'Second order line Curve', 'LineWidth', 1.5,
'MarkerSize', 10)
hold on
plot(n,Zi,'b-p', 'DisplayName', 'First order line Curve', 'LineWidth', 1.5,
'MarkerSize', 10)
hold on
plot(n,x1,'bo','DisplayName', 'Xi1 (Blue Points)', 'LineWidth', 1.5)
hold on
plot(n,x2,'ro','DisplayName', 'Xi2 (Red Points)', 'LineWidth', 1.5)
hold on
plot(n,y,'g*','DisplayName', 'Yi (Star Points)', 'LineWidth', 1.5)
title('Comparision of First order and Second order regression with Multivariable');
xlabel('Data Point Index (X)');
ylabel('Values (Y)');
grid on;
legend('show', 'Location', 'NorthWest');
hold off

*Chart is attached separately.

Result:
A Second order fit line equation is :
Yi=2.218019 + 2.493136 X1+ 2.004790 X2+0.141261 X1^2 -0.727464 X2^2
A first order best fit line equation is :
Zi=-3.014407 + 4.967827 X1 -1.782482 X2
Coefficient of Correlation is : 0.918558
Coefficient of Correlation is : 0.963199

Conclusion:
The second-order fit model successfully captures the underlying non-linear trends
in the data, resulting in a much stronger relationship between the independent
variables X1 and X2 and the dependent variable (Yi).
The second-order model has a Coefficient of Correlation R of 0.963199. This value
is substantially higher and closer to the ideal value of 1.0 when compared to the
first-order model's R value of 0.918558.
 Tutorial 10: Multivariable Second Order Regression

Problem: Find the best fit line curve and coefficient of correlation for following
data for multivariable second order regression.
X1 0 5 4 1 4
X2 0 1 2 3 6
Yi 1 7 12 8 3

Description: Multivariable Second-Order Regression (also known as a Quadratic

Model or a specific type of Multiple Linear Regression). It's a foundational
technique used in engineering and statistics to model curved or non-linear
relationships between a single output (response) variable and two or more input
variables (factors).

Flow Chart:

Start

Initialize the matrix X1, X2

and Y

Arrange in the form of

Y=X*A

A=Y/X
A=(XT*Y)/XT*X)

Find A,St,Sr and R

Display
Results
Main Program:
%% Multivariable Second Order Regression and Correlation Calculation
clc;
clear all;
Yi= [1; 7; 12; 8; 3; 6; 10; 5];
X1=[0; 4; 5; 1; 4; 2; 6; 3];
X2=[0; 1; 8; 15; 6; 4; 10; 7];
n=length(X1);
ones_column = ones(n, 1);
X_X1_sq = X1.^2;
X_X2_sq = X2.^2;
X = [ones_column, X1, X2, X_X1_sq, X_X2_sq];
B = X' * Yi;
C = X' * X;
A = C \ B;
y_bar = mean(Yi);
St = sum((Yi - y_bar).^2);
Y_pred = X * A;
Sr = sum((Yi - Y_pred).^2);
R = sqrt((St-Sr)/St);
disp(' REGRESSION ANALYSIS RESULTS');
disp('Regression Coefficients (A):');
labels = {'Intercept (A(1))', 'Coefficient X1 (A(2))', 'Coefficient X2 (A(3))', ...
'Coefficient X1^2 (A(4))', 'Coefficient X2^2 (A(5))', 'Coefficient X1*X2
(A(6))'};
for i = 1:length(A)
 fprintf('%-25s: %8.4f\n', labels{i}, A(i));
end
disp('-------------------------------------');
fprintf('Total Sum of Squares (St) : %8.4f\n', St);
fprintf('Residual Sum of Squares (Sr) : %8.4f\n', Sr);
fprintf('Coefficient of Correlation (R) : %8.4f\n', R);

Result:
REGRESSION ANALYSIS RESULTS

Regression Coefficients (A):

Intercept (A(1)) : 1.5086

Coefficient X1 (A(2)) : 1.3935

Coefficient X2 (A(3)) : -0.2607

Coefficient X1^2 (A(4)) : -0.0115

Coefficient X2^2 (A(5)) : 0.0400

-------------------------------------

Total Sum of Squares (St) : 90.0000

Residual Sum of Squares (Sr) : 31.6201

Coefficient of Correlation (R) : 0.8054

Conclusion:
Based on the calculated coefficients, the fitted regression equation is:
Y = 1.5086 + 1.3935 X1 - 0.2607 X2 - 0.0115 X12 + 0.0400 X22
Coefficient of Correlation (R): 0.8054 this value indicates a strong positive
correlation between the combined set of predictor variables (the linear and
quadratic terms of X1 and X2 and the response Y. The second-order model provides
a strong fit to the data, successfully explaining a significant majority of the
observed variability in the response.
 Tutorial 11: Gauss Newton Method

Problem: Problem: Find the best-fit non-linear curve and the final parameters
(a0,a1) for the following data using the Gauss-Newton method, fitting the model:
y = a0(1 – e(-a1*x))
X 0.25 0.75 1.25 1.75 2.25
Y 0.28 0.57 0.68 0.74 0.79

Description: The Gauss-Newton Method is an iterative technique used to find the

best-fit parameters for a non-linear model. Unlike linear regression, which has a
single closed-form solution, this method repeatedly refines an initial guess of the
parameters (a) until the Sum of Squared Errors (SSE) between the data (Y) and the
model (f) is minimized.

Flow Chart:
Start

Initialize the matrix X and Y and number of iteration

Arrange in the form of

D=Z*A

A=D/Z, A=(ZT*D)/ZT*Z)

Find A

No
If iteration is
less than n

Yes

Display
Results
Main Program:
% Gauss Newton Method
clc;
clear all;
% --- Initialize Data ---
X = [0.25; 0.75; 1.25; 1.75; 2.25];
Y = [0.28; 0.57; 0.68; 0.74; 0.79];
N = length(X);
% --- Initial Guess ---
a = [2; 2]; % a(1) = a0 (Asymptote), a(2) = a1 (Rate)
num_iterations = 5;

disp('GAUSS-NEWTON METHOD FOR NONLINEAR REGRESSION');

for iter = 1:num_iterations

f = a(1) * (1 - exp(-a(2) * X));

D = Y - f;
df = 1 - exp(-a(2) * X);
df2 = a(1) * X .* exp(-a(2) * X);

Z = [df, df2];
Z_T_Z = Z' * Z;
Z_T_r = Z' * D;
da = Z_T_Z \ Z_T_r;
 a = a + da;
SSE = sum(D.^2);

disp(['Iteration ', num2str(iter), ':']);

disp([' a = [', num2str(a(1), 4), '; ', num2str(a(2), 4), ']']);

end

% --- Final Fit Calculation (for R-squared) ---

Y_pred = a(1) * (1 - exp(-a(2) * X));
Sr = sum((Y - Y_pred).^2);

Y_bar = mean(Y);
St = sum((Y - Y_bar).^2);

R = sqrt((St - Sr) / St);

disp('FINAL RESULTS AFTER 5 ITERATIONS:');

disp(['a0 = ', num2str(a(1), 6)]);
disp(['a1 = ', num2str(a(2), 6)]);
disp('----------------------------------');
fprintf('Sum of Squares of Residuals (Sr): %8.6f\n', Sr);
fprintf('Coefficient of Correlation (R): %8.6f\n', R);
figure;
hold on;
plot(X, Y, 'o', 'MarkerFaceColor', 'b', 'MarkerSize', 8);

% Plot the fitted function

x_fit = linspace(min(X), max(X) + 1, 100);
y_fit = a(1) * (1 - exp(-a(2) * x_fit));

plot(x_fit, y_fit, 'r-', 'LineWidth', 2);

xlabel('x');
ylabel('y');
title('Gauss-Newton Nonlinear Regression Fit: y = a_0(1 - e^{-a_1 x})');
legend('Original Data', 'Best Fit Curve', 'Location', 'SouthEast');
grid on;
hold off;

Result:
GAUSS-NEWTON METHOD FOR NONLINEAR REGRESSION
Iteration 1:
a = [0.7863; 1.873]
Iteration 2:
a = [0.7891; 1.68]
Iteration 3:
a = [0.7918; 1.675]
Iteration 4:
a = [0.7919; 1.675]
Iteration 5:
a = [0.7919; 1.675]
FINAL RESULTS AFTER 5 ITERATIONS:
a0 = 0.791868
a1 = 1.67514
----------------------------------
Sum of Squares of Residuals (Sr): 0.000662
Coefficient of Correlation (R): 0.997989

Conclusion:
The Gauss-Newton method successfully converged over 5 iterations, using the
initial guess a = [2; 2], to find the optimal parameters for the non-linear growth
model.
The best-fit curve is described by the equation:
Y = 0.7919 (1 - e(-1.6751 X))
Here, coefficient of correlation is 0.997989 which states that the Gauss-Newton
method yielded an exceptionally accurate and stable non-linear model,
demonstrating excellent predictive capability for this dataset.
 Tutorial 12 : Romberg Extrapolation

Problem: Solve the following equation by using the romberg extrapolation

method. f(x)= 0.2 + 25x - 200x2 + 675x3 - 900x4 + 400x5 for 0 to 0.8.

Description: While basic methods like the Trapezoidal Rule are simple, they often
require a very high number of segments to achieve high precision. The goal of this
program is to implement a more efficient algorithm that achieves higher-order
accuracy with fewer function evaluations by combining the Trapezoidal Rule with
Romberg Integration.

Flow chart: -

START

Create an n×n matrix R

and calculate the initial
step size h.

For i = 1 to n
Calculate 2^{i-1} segments.
If i=1, calculate the basic trapezoid.
If i>1, calculate new points and update
R(i, 1) based on R(i-1, 1).

Apply Romberg Extrapolation to refine

previous estimates.

STOP
Main Program:
function [I, R] = romberg_integration(f, a, b, n)
clc;

if nargin == 0
romberg_test();
return;
end

if n < 1
error('N (number of levels) must be an integer greater than or equal to 1.');
end

R = zeros(n, n);
h = b - a;

for i = 1:n
num_segments = 2^(i-1);
h_i = (b - a) / num_segments;

if i == 1
R(i, 1) = h_i * (f(a) + f(b)) / 2;
else
sum_new_points = 0;
 for j = 1:2:(num_segments - 1)
x = a + j * h_i;
sum_new_points = sum_new_points + f(x);
end

R(i, 1) = 0.5 * R(i-1, 1) + h_i * sum_new_points;

end
end

for j = 2:n
for i = j:n
power_of_4 = 4^(j-1);

R(i, j) = (power_of_4 * R(i, j-1) - R(i-1, j-1)) / (power_of_4 - 1);

end
end

I = R(n, n);
disp('ROMBERG INTEGRATION RESULTS');
disp(['Function: ', func2str(f), ', Limits: [', num2str(a), ', ', num2str(b), '], Levels: ',
num2str(n)]);
disp(' ');
disp('Complete Romberg Table (R):');
disp(R);
disp(' ');
disp(['Final Integral Estimate (R(', num2str(n), ',', num2str(n), ')): ', num2str(I,
10)]);

end
function romberg_test()
disp('*** Running Example Demonstration ***');

f_example = @(x) 0.2 + 25*x - 200*x.^2 + 675*x.^3 - 900*x.^4 + 400*x.^5;

a_example = 0;
b_example = 0.8;
n_example = 7;

[I_final, ~] = romberg_integration(f_example, a_example, b_example,

n_example);
analytical_result = 1.640533333333333;
disp(' ');
disp('--- Analytical Comparison ---');
disp(['Analytical Result: ', num2str(analytical_result, 10)]);
disp(['Numerical Error: ', num2str(abs(I_final - analytical_result), 10)]);
disp('--- Example Execution Complete ---');
end

Result:
ROMBERG INTEGRATION RESULTS

Function: @(x)0.2+25*x-200*x.^2+675*x.^3-900*x.^4+400*x.^5, Limits: [0, 0.8], Levels: 7

Complete Romberg Table (R):

 0.1728 0 0 0 0 0 0

1.0688 1.3675 0 0 0 0 0

1.4848 1.6235 1.6405 0 0 0 0

1.6008 1.6395 1.6405 1.6405 0 0 0

1.6306 1.6405 1.6405 1.6405 1.6405 0 0

1.6380 1.6405 1.6405 1.6405 1.6405 1.6405 0

1.6399 1.6405 1.6405 1.6405 1.6405 1.6405 1.6405

Final Integral Estimate (R(7,7)): 1.640533333

--- Analytical Comparison ---

Analytical Result: 1.640533333

Numerical Error: 1.33226763e-15

Conclusion: The method converged rapidly and essentially provided the exact
value. The Romberg Integration method successfully solved the equation, yielding
a result that is nearly exact due to the function being a polynomial. We got the
result in just n=3 (2-Segment).The convergence to the analytical answer is
achieved quickly in the table, demonstrating the high efficiency of the method.