Lagrange interpolation and Newton's interpolation - MATLAB implementation

Lagrange interpolation and Newton's interpolation

Special paper of undergraduate experiment report

拉格朗日插值，牛顿差值，MATLAB

Experimental purpose and requirement

Experimental Purpose: The purpose of this experiment is to compare the accuracy and efficiency of Lagrange interpolation and Newton's interpolation methods for approximating a function based on a set of discrete data points. The two interpolation methods will be implemented and their results will be compared based on the error between the approximated values and the actual function values.

Experimental Requirements: To conduct this experiment, the following requirements must be fulfilled:

  1. Data Collection: A set of discrete data points will be collected to approximate a given function. Here we use data on the textbook and some randomised data.
  2. Lagrange Interpolation: The Lagrange interpolation method will be implemented using the collected data points. The implementation should be able to take in a set of data points and approximate the function at any given point.
  3. Newton's Interpolation: The Newton's interpolation method will also be implemented using the collected data points. The implementation should be able to take in a set of data points and approximate the function at any given point.
  4. Error Calculation: The error between the approximated values and the actual function values will be calculated using a suitable error metric such as the Mean Squared Error or Absolute Error.
  5. Comparison of Results: The results obtained from both the interpolation methods will be compared based on the error between the approximated values and the actual function values. The accuracy and efficiency of the two methods will be evaluated based on this comparison.

Experimental principle and main content

Equipment, Toolbox

Only Symbolic Math Toolbox is used here.

It can be done without using Symbolic Math Toolbox, yielding a numerical value at the end, but this is not implemented here. In this report, we first calculate the interpolated polynomial symbolically, then compute their numerical value.

Core formulae: Lagrange interpolation

The Lagrange interpolated polynomial is given by

\[p=\sum_{k=0}^n f\left(x_k\right) l_k \]

where

\[l_k(x)=\prod_{\substack{j=0 \\ j \neq k}}^n\left(x-x_j\right) /\left(x_k-x_j\right), \quad a \leqslant x \leqslant b. \]

Core formulae: Newton's interpolation

The main idea of Newton's interpolation is

\[p_{k+1}(x)=p_k(x)+\left\{\prod_{j=0}^k\left(x-x_j\right)\right\} f\left[x_0, x_1, \ldots, x_{k+1}\right], \quad a \leqslant x \leqslant b \]

which is the same as

\[\begin{aligned} p_n(x)= & f\left(x_0\right)+\left(x-x_0\right) f\left[x_0, x_1\right]+\left(x-x_0\right)\left(x-x_1\right) f\left[x_0, x_1, x_2\right] \\ & +\ldots+\left\{\prod_{j=0}^{n-1}\left(x-x_j\right)\right\} f\left[x_0, x_1, \ldots, x_n\right], \quad a \leqslant x \leqslant b, \end{aligned} \]

And the dividend difference term is recursively calculate by

\[f\left[x_j, x_{j+1}, \ldots, x_{j+k+1}\right]=\frac{f\left[x_{j+1}, \ldots, x_{j+k+1}\right]-f\left[x_j, \ldots, x_{j+k}\right]}{\left(x_{j+k+1}-x_j\right)} \]

Before designed our code, we can do some simple calculation to check whether our thought is correct.

The recurrence relation \((5.14)\) was used to calculate the divided differences of the function

\[f(x)=10 \mathrm{e}^{-3 x}, \quad 0 \leqslant x \leqslant 2, \]

tabulated on the point set \(\{1.60,1.63,1.70,1.76,1.80\}\). The results are shown in Tables 5.1 and 5.2. All data and all calculated numbers were rounded to a fixed precision before they were recorded and used for subsequent calculation. The difference between the tables is that in Table

By the data in table 5.2, we have that

\[p = 0.0823 \\- 0.23633 (-1.6 + x) + 0.329 (-1.63 + x) (-1.6 + x) - \\ 0.32887 (-1.7 + x) (-1.63 + x) (-1.6 + x) + \\ 0.5008 (-1.76 + x) (-1.7 + x) (-1.63 + x) (-1.6 + x)\\ = 0.5008 x^4-3.67922 x^3+10.3516 x^2-13.3214 x+6.68435 \]

By considering rounding,

clear;
% range_x = [-10, 10];
% sz_vec = [4, 1];
% rng(1);
% x_vec = randperm(range_x(2), sz_vec(1)).';
% y_vec = randperm(range_x(2), sz_vec(1)).';
x_vec = [1.60, 1.63, 1.70, 1.76, 1.80];
y_vec = 10.* exp(-3.*x_vec);
preci = 5;
n = length(x_vec);
F = zeros([n, n]);
for i = 1:n
    F(i,i) = y_vec(i);
    F(i,i) = round(F(i,i), preci);
end

% calculate the divided difference
for delta = 1:(n-1)
    for i = 1:n
    
    
        if i + delta > n 
            break;
        end
        j = i + delta;
        F(i, j) = (F(i, j-1) - F(i+1, j)) ./ (x_vec(i) - x_vec(j));
        F(i, j) = round(F(i, j), preci);
    end
end


syms x
p = 0;
cur_prod = 1;
for i = 1:n
    if i >= 2
        cur_prod = cur_prod * (x - x_vec(i - 1));
    end
    p = p + cur_prod * F(1, i);
    % disp(p);
    
end

p = simplify(p);
disp(latex(p))
p = matlabFunction(p);

min_num = min(x_vec);
max_num = max(x_vec);
total_value = 1000;
xi = linspace(min_num, max_num, total_value);
plot(x_vec, y_vec, 'o', xi, p(xi))
legend('Data points', 'Lagrange interpolation')

\[p = \frac{10017\,x^4}{20000}-\frac{7359133\,x^3}{2000000}+\frac{1035253783\,x^2}{100000000}-\frac{1665300473\,x}{125000000}+\frac{1044497349}{156250000} \]

(which is the same as)

p = x.*(-1.3322403784e+1)+x.^2.*1.035253783e+1-x.^3.*3.6795665+x.^4.*5.0085e-1+6.6847830336

This result is nearly the same as in the book.

1_eg

Experimental result and analysis

Test of interpolation

Random different points

test1 
>> all_in_one
(31*x^9)/20160 - (127*x^8)/1680 + (2293*x^7)/1440 - (2707*x^6)/144 + (393197*x^5)/2880 - (226591*x^4)/360 + (2304703*x^3)/1260 - (201361*x^2)/63 + (11981*x)/4 - 1111
symbolic function inputs: x
 
Elapsed time is 0.406410 seconds.
2*x + ((x - 8)*(x - 9))/2 - (2*(x - 7)*(x - 8)*(x - 9))/3 - (29*(x - 7)*(x - 8)*(x - 9)*(x - 10))/280 - ((x - 3)*(x - 7)*(x - 8)*(x - 9)*(x - 10))/210 + (239*(x - 3)*(x - 4)*(x - 7)*(x - 8)*(x - 9)*(x - 10))/5040 + (103*(x - 3)*(x - 4)*(x - 6)*(x - 7)*(x - 8)*(x - 9)*(x - 10))/10080 + (13*(x - 2)*(x - 3)*(x - 4)*(x - 6)*(x - 7)*(x - 8)*(x - 9)*(x - 10))/10080 + (31*(x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 6)*(x - 7)*(x - 8)*(x - 9)*(x - 10))/20160 - 13
symbolic function inputs: x
 
Elapsed time is 0.198716 seconds.

It can be found that the oscillation is more obvious at the border.

Some real functions

sin(x)
>> all_in_one
(sin(2)*(x - 1)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 7)*(x - 8)*(x - 9)*(x - 10))/40320 - (sin(1)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 7)*(x - 8)*(x - 9)*(x - 10))/362880 - (sin(3)*(x/2 - 1/2)*(x - 2)*(x - 4)*(x - 5)*(x - 6)*(x - 7)*(x - 8)*(x - 9)*(x - 10))/5040 + (sin(4)*(x/3 - 1/3)*(x - 2)*(x - 3)*(x - 5)*(x - 6)*(x - 7)*(x - 8)*(x - 9)*(x - 10))/1440 - (sin(5)*(x/4 - 1/4)*(x - 2)*(x - 3)*(x - 4)*(x - 6)*(x - 7)*(x - 8)*(x - 9)*(x - 10))/720 + (sin(6)*(x/5 - 1/5)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 7)*(x - 8)*(x - 9)*(x - 10))/576 - (sin(7)*(x/6 - 1/6)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 8)*(x - 9)*(x - 10))/720 + (sin(8)*(x/7 - 1/7)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 7)*(x - 9)*(x - 10))/1440 - (sin(9)*(x/8 - 1/8)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 7)*(x - 8)*(x - 10))/5040 + (sin(10)*(x/9 - 1/9)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 7)*(x - 8)*(x - 9))/40320
symbolic function inputs: x
 
Elapsed time is 0.939289 seconds.
(76365784749291*x)/1125899906842624 - (7530053351810631*(x - 1)*(x - 2))/18014398509481984 + (8481884696636481*(x - 1)*(x - 2)*(x - 3))/72057594037927936 + (1432589426173745*(x - 1)*(x - 2)*(x - 3)*(x - 4))/288230376151711744 - (7292257428457373*(x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5))/1152921504606846976 + (941795826732659*(x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6))/1152921504606846976 + (4599379538442441*(x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 7))/147573952589676412928 - (2507741590142589*(x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 7)*(x - 8))/147573952589676412928 + (1579577787439489*(x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 7)*(x - 8)*(x - 9))/1180591620717411303424 + 3484185274626763/4503599627370496
symbolic function inputs: x
 
Elapsed time is 0.217115 seconds.
test2_real
floor(log(x)), range from 1 to 10
>> all_in_one
(x*(5*x^8 - 225*x^7 + 4260*x^6 - 44100*x^5 + 271131*x^4 - 1003275*x^3 + 2145460*x^2 - 2354400*x + 981144))/90720
symbolic function inputs: x
 
Elapsed time is 0.262802 seconds.
max(|f - g|) = 1.4724
var(|f - g|) = 0.11019
Elapsed time is 0.054497 seconds.
((x - 1)*(x - 2))/2 - ((x - 1)*(x - 2)*(x - 3))/3 + ((x - 1)*(x - 2)*(x - 3)*(x - 4))/8 - ((x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5))/30 + ((x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6))/144 - ((x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 7))/1008 + ((x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 7)*(x - 8)*(x - 9))/18144
symbolic function inputs: x
 
max(|f - g|) = 1.4724
var(|f - g|) = 0.11019

In this case, no significant error appears.

floorlog
floor(log(x)), range from 1 to 100
testlarge

Lagrange interpolation:

Elapsed time is 19.243801 seconds.
max(|f - g|) = 2.791190982682177e+25
var(|f - g|) = 4.164987473892946e+48

Newton's interpolation:

Elapsed time is 0.226222 seconds.
max(|f - g|) = 8.322704101592106e+26
var(|f - g|) = 3.989048826417586e+51

In this case, Lagrange interpolation is slower but have a subtlety higher precision.

1/(1+x^2)
>> all_in_one
- x^10/44200 + (7*x^8)/5525 - (83*x^6)/3400 + (2181*x^4)/11050 - (149*x^2)/221 + 1
symbolic function inputs: x
 
Elapsed time is 0.328900 seconds.
max(|f - g|) = 1.9156
var(|f - g|) = 0.25155
Elapsed time is 0.032546 seconds.
(9*x)/442 + (23*(x + 4)*(x + 5))/2210 + (7*(x + 3)*(x + 4)*(x + 5))/1105 + (19*(x + 2)*(x + 3)*(x + 4)*(x + 5))/4420 - (9*(x + 1)*(x + 2)*(x + 3)*(x + 4)*(x + 5))/4420 - (x*(x + 1)*(x + 2)*(x + 3)*(x + 4)*(x + 5))/884 + (6*x*(x - 1)*(x + 1)*(x + 2)*(x + 3)*(x + 4)*(x + 5))/5525 - (19*x*(x - 1)*(x + 1)*(x - 2)*(x + 2)*(x + 3)*(x + 4)*(x + 5))/44200 + (x*(x - 1)*(x + 1)*(x - 2)*(x + 2)*(x - 3)*(x + 3)*(x + 4)*(x + 5))/8840 - (x*(x - 1)*(x + 1)*(x - 2)*(x + 2)*(x - 3)*(x + 3)*(x - 4)*(x + 4)*(x + 5))/44200 + 31/221
symbolic function inputs: x
 
max(|f - g|) = 1.9156
var(|f - g|) = 0.25155
cheb

Conclusion

In this experiment, we have implemented two methods for polynomial interpolation: Lagrange interpolation and Newton's interpolation. Both methods are useful for approximating a function when only a limited number of data points are available.

Lagrange interpolation involves constructing a polynomial that passes through all of the given data points. This polynomial is determined by using a set of basis polynomials, which are then weighted by scalar coefficients to produce the final polynomial. One advantage of Lagrange interpolation is that the resulting polynomial is unique, and therefore it can be used to perform a more accurate predict value of the function at any point within the range of the data points.

Newton's interpolation, on the other hand, involves constructing a polynomial using divided differences. This method allows for the efficient calculation of the polynomial coefficients, and therefore it can be used to interpolate functions even when the number of data points is large. Additionally, Newton's interpolation can be extended to handle unevenly spaced data points, which makes it a more versatile method than Lagrange interpolation in some cases.

Overall, both Lagrange interpolation and Newton's interpolation are powerful methods for approximating functions, and the choice between the two will depend on the specific requirements of the problem at hand.

Code


% initialization
clear;
preci = 10;
syms x

f(x) = 1/(1+x^2);
f_range = [-5, 5];
[x_vec, y_vec] = init(2, f, f_range);
tiledlayout(2,1);
total_value = 1000;

% lagrange interpolation
tic
p(x) = lagrange_interpolation(x_vec, y_vec);
disp(p);
toc
nexttile
show_graph(p, x_vec, y_vec, 1000, "lagrange interpolation");
error_analysis(f, p, f_range, total_value);

% newton's interpolation
tic 
F = newtons_interpolation_init(x_vec, y_vec);
p(x) = newtons_interpolation(x_vec, F);
toc
disp(p);
nexttile
show_graph(p, x_vec, y_vec, 1000, "newtons interpolation");
error_analysis(f, p, f_range, total_value);

function [x_vec, y_vec] = init(init_type, f, f_range)
    if init_type == 1
        range_x = [0, 10];
        sz_vec = [10, 1];
        % randperm initialization
        x_vec = randperm(range_x(2), sz_vec(1)).';
        y_vec = randperm(range_x(2), sz_vec(1)).';
    end
    if init_type == 2
        % real function initialization
        f_min = min(f_range);
        f_max = max(f_range);
        x_vec = f_min:f_max;
        x_vec = x_vec.';
        y_vec = f(x_vec);
    end 
    if init_type == 3
        % example in textbook
        x_vec = [1.60, 1.63, 1.70, 1.76, 1.80];
        y_vec = 10 .* exp(-3 .* x_vec);
    end
end

function p = lagrange_interpolation(x_vec, y_vec)
    % Define the Lagrange basis polynomials
    syms x;
    L = sym('L', [1,length(x_vec)]);
    
    n = length(x_vec);
    for i = 1:n
        L(i) = 1;
        for j = 1:n
            if j ~= i
                L(i) = L(i) * (x - x_vec(j)) / (x_vec(i) - x_vec(j));
            end
        end
    end
    p = simplify(sum(y_vec.' .* L));
end

function F = newtons_interpolation_init(x_vec, y_vec)
    n = length(y_vec);
    F = zeros([n, n]);
    
    % start to count time 
    tic
    
    % calculate the divided difference
    for i = 1:n
        F(i,i) = y_vec(i);
        % F(i,i) = round(F(i,i), preci);
    end
    
    for delta = 1:(n-1)
        for i = 1:n
            if i + delta > n 
                break;
            end
            j = i + delta;
            F(i, j) = (F(i, j-1) - F(i+1, j)) ./ (x_vec(i) - x_vec(j));
            % F(i, j) = round(F(i, j), preci);
        end
    end
end

function p = newtons_interpolation(x_vec, F)
    % calculate the interpolation polynomial
    n = length(x_vec);
    syms x
    p = 0;
    cur_prod = 1;
    for i = 1:n
        if i >= 2
            cur_prod = cur_prod * (x - x_vec(i - 1));
        end
        p = p + cur_prod * F(1, i);
        % disp(p);
    end
end

function show_graph(p, x_vec, y_vec, total_value, str)
    min_num = min(x_vec);
    max_num = max(x_vec);
    xi = linspace(min_num, max_num, total_value);
    plot(x_vec, y_vec, 'o', xi, p(xi));
    legend('data points', str);
end

function error_analysis(f, g, f_range, total_value)
    f_min = min(f_range);
    f_max = max(f_range);
    xi = linspace(f_min, f_max, total_value);
    ei = double(abs(f(xi) - g(xi)));
    max_dev = max(ei);
    var_dev = var(ei);
    disp("max(|f - g|) = " + num2str(max_dev));
    disp("var(|f - g|) = " + num2str(var_dev));

end
posted @ 2023-04-13 14:07  miyasaka  阅读(72)  评论(0)    收藏  举报