Numerical Methods Revision

🧠 Part 1: Core Concepts & Formulas (Quick Review)

🔹 Root Finding (TP3)
Method Formula / Update Rule Stopping Convergence
Criterion
Bisection c = , replace a or b
a+b

2
abs(b-a) < ε Linear, guaranteed if
based on sign of f (c) f (a)f (b) < 0

Newton x k+1 = x k −
f (x k )
′
f (x k )
abs(x_new - Quadratic (fast),
x_old) < ε needs f (x) ≠ 0
′

Fixed x k+1 = g(x k ) abs(x_new - Linear if g'(x)<1 near

Point x_old) < ε the root

Theoretical Bisection Iterations:

b−a
ln ( )
ε
N max ≥
ln(2)

🔹 Interpolation (TP4 Part A & B)

Interpolation: Polynomial passes exactly through all given nodes.
Lagrange: P , where L
n x−x j
n (x) = ∑ i=0 y i L i (x) i (x) = ∏ j≠i
x i −x j

Newton: P n (x) = f [x 0 ] + f [x 0 , x 1 ](x − x 0 ) + f [x 0 , x 1 , x 2 ](x − x 0 )(x − x 1 ) + …

Divided Differences: f [x
f [x i+1 ,…,x i+k ]−f [x i ,…,x i+k−1 ]
i, … , x i+k ] =
x i+k −x i

Lagrange and Newton produce the same unique polynomial (just different
representations).

🔹 Approximation (TP4 Part C)

Least Squares: Finds coefficients that minimize ∑(y i − y
^i )
2
. Does not necessarily pass
through points.
Linear LS: a = n ∑ x i y i −∑ x i ∑ y i
2
n ∑ x −(∑ x i )
2
,b= ∑ y i −a ∑ x i

n
i

Polynomial LS: Solve normal equations: (A T

A)c = A
T
y

💻 Part 2: Code Templates

How to use: Copy each block into a separate .m file or paste directly into the command
window. Change x , y , or tol as needed for test variations.
📘 TP3 Exercise 1: Bisection Method
% TP3_Ex1_Bisection.m
clear; clc;

f = @(x) x.^3 - x - 2;
a = 1; b = 2;
tol = 1e-6;

% 1. Verify sign change

if f(a)*f(b) > 0
error('No root guaranteed in [a,b]. Check f(a)*f(b) < 0.');
end

% 2. Theoretical max iterations

N_theory = ceil( log((b-a)/tol) / log(2) );
fprintf('Theoretical max iterations: %d\n', N_theory);

% 3-6. Bisection loop

iter = 0;
x_history = [];
while abs(b-a) > tol
c = (a+b)/2;
iter = iter + 1;
x_history = [x_history; c];

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

fprintf('Bisection iterations: %d\n', iter);

fprintf('Approximated root: %.8f\n', c);

% 7. Convergence plot
figure;
plot(1:iter, x_history, 'b-o', 'LineWidth', 1.5, 'MarkerSize', 5);
grid on;
xlabel('Iteration');
ylabel('Approximation x_k');
title('Bisection Method Convergence');

📘 TP3 Exercise 2: Newton's Method

 % TP3_Ex2_Newton.m
clear; clc;

f = @(x) x.^3 - x - 2;
df = @(x) 3*x.^2 - 1; % Derivative

x0 = 1.5;
tol = 1e-6;

iter = 0;
x_old = x0;
x_history = [];

while true
iter = iter + 1;
x_new = x_old - f(x_old)/df(x_old);
x_history = [x_history; x_new];

if abs(x_new - x_old) < tol

break;
end
x_old = x_new;
end

fprintf('Newton iterations: %d\n', iter);

fprintf('Approximated root: %.8f\n', x_new);

figure;
plot(1:iter, x_history, 'r-s', 'LineWidth', 1.5, 'MarkerSize', 5);
grid on;
xlabel('Iteration');
ylabel('Approximation x_k');
title('Newton Method Convergence');

📘 TP3 Exercise 3: Fixed Point Method

% TP3_Ex3_FixedPoint.m
clear; clc;

g = @(x) (x + 2).^(1/3);
x0 = 1;
tol = 1e-6;

iter = 0;
x_old = x0;
x_history = [];

while true
 iter = iter + 1;
x_new = g(x_old);
x_history = [x_history; x_new];

if abs(x_new - x_old) < tol

break;
end
x_old = x_new;
end

fprintf('Fixed Point iterations: %d\n', iter);

fprintf('Approximated root: %.8f\n', x_new);

figure;
plot(1:iter, x_history, 'g-^', 'LineWidth', 1.5, 'MarkerSize', 5);
grid on;
xlabel('Iteration');
ylabel('Approximation x_k');
title('Fixed Point Method Convergence');

📘 TP3 Exercise 4: Comparison Plot

% TP3_Ex4_Comparison.m
% Run Ex1, Ex2, Ex3 first to get x_history_bis, x_history_new,
x_history_fix
% For standalone use, re-define them here:
% (Omitted for brevity; just call plot commands with saved histories)
figure;
plot(1:length(x_history_bis), x_history_bis, 'b-o', 'DisplayName',
'Bisection');
hold on;
plot(1:length(x_history_new), x_history_new, 'r-s', 'DisplayName',
'Newton');
plot(1:length(x_history_fix), x_history_fix, 'g-^', 'DisplayName', 'Fixed
Point');
grid on;
xlabel('Iteration');
ylabel('x_k');
title('Convergence Comparison');
legend('show');
hold off;

📘 TP4 Ex1: Lagrange Interpolation (Experimental Data)

% TP4_Ex1_Lagrange_Data.m
clear; clc;

x = [0, 1, 2, 3];
 y = [1, 3, 2, 5];
n = length(x) - 1;

% Function to evaluate Lagrange polynomial at any X

P = @(X) arrayfun(@(xx) sum( y .* arrayfun(@(i) prod( (xx - x([1:i-1,
i+1:n+1])) ./ (x(i) - x([1:i-1, i+1:n+1])) ), 1:n+1) ), X);

% 1. Plot points
figure;
plot(x, y, 'ko', 'MarkerFaceColor', 'k', 'MarkerSize', 8, 'DisplayName',
'Data');
hold on;

% 2-3. Plot polynomial

xx = linspace(min(x), max(x), 500);
yy = P(xx);
plot(xx, yy, 'b-', 'LineWidth', 2, 'DisplayName', 'Lagrange P_3(x)');

% 4. Estimate at x=1.5
x_est = 1.5;
y_est = P(x_est);
fprintf('P(%.1f) = %.4f\n', x_est, y_est);

% 5. Compare
plot(x_est, y_est, 'rv', 'MarkerSize', 10, 'MarkerFaceColor', 'r',
'DisplayName', sprintf('P(%.1f)', x_est));
grid on;
xlabel('x'); ylabel('y');
title('Lagrange Interpolation (Ex1)');
legend('show');
hold off;

📘 TP4 Ex2: Lagrange Interpolation of Known Function

% TP4_Ex2_Lagrange_Function.m
clear; clc;

f = @(x) log(1+x);
nodes = [0, 0.5, 1, 2];
y_nodes = f(nodes);
n = length(nodes)-1;

% Lagrange evaluator
LagrangeEval = @(X, xn, yn) arrayfun(@(xx) sum( yn .* arrayfun(@(i) prod(
(xx - xn([1:i-1, i+1:n+1])) ./ (xn(i) - xn([1:i-1, i+1:n+1])) ), 1:n+1) ),
X);

xx = linspace(min(nodes), max(nodes), 400);

yy_interp = LagrangeEval(xx, nodes, y_nodes);
 yy_exact = f(xx);

% Error
abs_error = abs(yy_exact - yy_interp);

figure;
subplot(2,1,1);
plot(nodes, y_nodes, 'ko', 'MarkerFaceColor', 'k'); hold on;
plot(xx, yy_exact, 'b-', 'LineWidth', 2, 'DisplayName', 'f(x)=ln(1+x)');
plot(xx, yy_interp, 'r--', 'LineWidth', 1.5, 'DisplayName', 'P_3(x)');
grid on; legend; title('Function vs Interpolation');

subplot(2,1,2);
plot(xx, abs_error, 'm-', 'LineWidth', 1.5);
grid on; xlabel('x'); ylabel('|f(x)-P(x)|');
title('Absolute Error');

📘 TP4 Ex3: Newton Interpolation (Divided Differences)

% TP4_Ex3_Newton_Data.m
clear; clc;

x = [1, 2, 4, 5];
y = [2, 3, 1, 4];
n = length(x);

% Build divided difference table

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

% Extract coefficients (first row)

c = DD(1,:);

% Horner evaluation for Newton form

NewtonEval = @(X) polyval(fliplr(c), X); % Note: polyval expects
descending powers, but Newton uses (x-x0)...
% Better explicit Horner for Newton:
NewtonEvalExact = @(X) arrayfun(@(xx) c(1) + c(2)*(xx-x(1)) + c(3)*(xx-
x(1))*(xx-x(2)) + c(4)*(xx-x(1))*(xx-x(2))*(xx-x(3)), X);

xx = linspace(min(x), max(x), 300);

yy = NewtonEvalExact(xx);

figure;
 plot(x, y, 'ko', 'MarkerFaceColor', 'k', 'DisplayName', 'Data'); hold on;
plot(xx, yy, 'b-', 'LineWidth', 2, 'DisplayName', 'Newton P(x)');
plot(3, NewtonEvalExact(3), 'rv', 'MarkerSize', 10, 'MarkerFaceColor',
'r', 'DisplayName', 'P(3)');
fprintf('P(3) = %.4f\n', NewtonEvalExact(3));
grid on; legend; xlabel('x'); ylabel('y'); title('Newton Interpolation
(Ex3)');

📘 TP4 Ex4: Newton Interpolation of Known Function

% TP4_Ex4_Newton_Function.m
clear; clc;

f = @(x) sqrt(x+1);
nodes = [0, 1, 2, 4];
y_nodes = f(nodes);
n = length(nodes);

DD = zeros(n, n);
DD(:,1) = y_nodes';
for j = 2:n
for i = 1:n-j+1
DD(i,j) = (DD(i+1,j-1) - DD(i,j-1)) / (nodes(i+j-1) - nodes(i));
end
end
c = DD(1,:);

% Newton evaluation
P_newton = @(X) arrayfun(@(xx) c(1) + c(2)*(xx-nodes(1)) + c(3)*(xx-
nodes(1))*(xx-nodes(2)) + c(4)*(xx-nodes(1))*(xx-nodes(2))*(xx-nodes(3)),
X);

xx = linspace(min(nodes), max(nodes), 400);

yy_exact = f(xx);
yy_interp = P_newton(xx);
err = abs(yy_exact - yy_interp);

figure;
plot(nodes, y_nodes, 'ko', 'MarkerFaceColor', 'k'); hold on;
plot(xx, yy_exact, 'b-', 'LineWidth', 2, 'DisplayName', 'f(x)=sqrt(x+1)');
plot(xx, yy_interp, 'r--', 'LineWidth', 1.5, 'DisplayName', 'Newton
P(x)');
grid on; legend; title('Newton Interpolation of f(x)');

figure;
plot(xx, err, 'm-', 'LineWidth', 1.5);
grid on; xlabel('x'); ylabel('Absolute Error'); title('Error
Distribution');
📘 TP4 Ex5: Linear Least Squares
% TP4_Ex5_LinearLS.m
clear; clc;

x = [1, 2, 3, 4, 5];
y = [2.2, 2.8, 3.6, 4.5, 5.1];
n = length(x);

% Coefficients
Sx = sum(x); Sy = sum(y);
Sxx = sum(x.^2); Sxy = sum(x.*y);

a = (n*Sxy - Sx*Sy) / (n*Sxx - Sx^2);

b = (Sy - a*Sx) / n;

fprintf('Linear fit: y = %.4f x + %.4f\n', a, b);

% Total squared error

y_pred = a*x + b;
SSE = sum((y - y_pred).^2);
fprintf('Total Squared Error: %.4f\n', SSE);

figure;
plot(x, y, 'ko', 'MarkerFaceColor', 'k', 'DisplayName', 'Data'); hold on;
xx = linspace(min(x), max(x), 100);
plot(xx, a*xx + b, 'r-', 'LineWidth', 2, 'DisplayName', 'Least Squares
Line');
grid on; legend; xlabel('x'); ylabel('y'); title('Linear Least Squares
Fit');

📘 TP4 Ex6: Quadratic Least Squares

% TP4_Ex6_QuadraticLS.m
clear; clc;

x = [0, 1, 2, 3, 4]';
y = [1.1, 2.0, 2.9, 3.7, 5.2]';
n = length(x);

% Build Vandermonde matrix A

A = [x.^2, x, ones(n,1)];

% Normal equations: (A'*A)*X = A'*y

X = (A'*A) \ (A'*y);
a = X(1); b = X(2); c = X(3);
fprintf('Quadratic fit: P(x) = %.4f x^2 + %.4f x + %.4f\n', a, b, c);
 % Plot
xx = linspace(min(x), max(x), 100);
yy_quad = a*xx.^2 + b*xx + c;

% Linear fit for comparison

Sx=sum(x); Sy=sum(y); Sxx=sum(x.^2); Sxy=sum(x.*y);
a_lin = (n*Sxy - Sx*Sy)/(n*Sxx - Sx^2);
b_lin = (Sy - a_lin*Sx)/n;
yy_lin = a_lin*xx + b_lin;

figure;
plot(x, y, 'ko', 'MarkerFaceColor', 'k', 'DisplayName', 'Data'); hold on;
plot(xx, yy_lin, 'b--', 'LineWidth', 1.5, 'DisplayName', 'Linear Fit');
plot(xx, yy_quad, 'r-', 'LineWidth', 2, 'DisplayName', 'Quadratic Fit');
grid on; legend; xlabel('x'); ylabel('y'); title('Linear vs Quadratic
Least Squares');

📝 Part 3: Reflection Answers

🔹 TP3 Reflections
Question Answer
Why does Bisection always It relies on the Intermediate Value Theorem. Each step
converge if f (a)f (b) < 0? halves the interval, guaranteeing the root stays inside.
Why is Newton faster? Quadratic convergence: error roughly squares each
step. Uses derivative to jump directly toward root.
When does Fixed Point fail? If the Derivative Condition is Violated.

🔹 TP4 Ex1–2 Reflections (Lagrange)

Question Answer
Why does polynomial pass By construction, L i
(x j ) = δ ij . So P (x j
) = yj .
exactly through all points?
What happens with many Runge's phenomenon: high-degree polynomials oscillate
points? wildly near edges. Condition number grows.
Interpolation vs Interpolation: exact at nodes, global polynomial.
Approximation? Approximation: best fit (e.g., least squares), minimizes
error, doesn't force exact matches.

🔹 TP4 Ex3–4 Reflections (Newton)

Question Answer
Main advantage of Easy to add new nodes without recomputing everything.
Newton's form? Only one new divided difference needed.
What happens if a new Append one row/column to divided difference table. Old
node is added? coefficients stay unchanged.
Why practical numerically? Builds incrementally, stable for adding points, easy to
evaluate via Horner-like scheme.

🔹 TP4 Ex5–6 Reflections (Least Squares)

Question Answer
Why doesn't LS line pass It minimizes sum of squared residuals. Exact match is
through all points? usually impossible with noisy/overdetermined data.
Difference interpolation vs Interpolation forces P (x ) = y . LS minimizes i i

LS? ∑(y − P (x )) , trades exactness for smoothness/stability.

i i
2

What does total squared Measure of overall deviation between model and data.
error represent? Lower = better fit.
Why quadratic fits better? Captures curvature. Linear assumes constant slope;
quadratic allows acceleration/change in slope.
What if degree is too high? Overfitting: fits noise, oscillates, poor prediction outside
data.

🗂️ Part 4: Cheat Sheet & Documentation

📐 Essential Formulas
Method Formula
Bisection step c =
a+b

Bisection iterations N ≥
ln((b−a)/ε)

ln 2

Newton update x k+1 = x k −

f (x k )
′
f (x k )

Fixed Point x k+1 = g(x k )

Lagrange basis L i (x) = ∏ j≠i

x−x j

x i −x j

Divided Difference f [x i , … , x i+k ] =

f [x i+1 ,…,x i+k ]−f [x i ,…,x i+k−1 ]

x i+k −x i

Linear LS slope a =
n ∑ xy−∑ x ∑ y
2 2
n ∑ x −(∑ x)

Linear LS intercept b =
∑ y−a ∑ x

n
Method Formula
Polynomial LS X = (A
T
A)
−1
A
T
y