La méthode d'interpolation de Lagrange

- Je utiliser la convolution et pour les boucles (trop pour les boucles) pour le calcul de l'interpolation à l'aide de
Lagrange's method , voici le code principal :

function[p] = lagrange_interpolation(X,Y)
L = zeros(n);
p = zeros(1,n);


% computing L matrice, so that each row i holds the polynom L_i
% Now we compute li(x) for i=0....n  ,and we build the polynomial 

for k=1:n
    multiplier = 1;
    outputConv = ones(1,1);
    for index = 1:n
        if(index ~= k && X(index) ~= X(k))
            outputConv = conv(outputConv,[1,-X(index)]);
            multiplier = multiplier * ((X(k) - X(index))^-1);
        end
    end
    polynimialSize = length(outputConv);
    for index = 1:polynimialSize
        L(k,n - index + 1) = outputConv(polynimialSize - index + 1);
    end
    L(k,:) = multiplier .* L(k,:);
end

% continues 


end

Ceux qui sont trop pour les boucles pour le calcul de la l_i(x) (ce qui est fait avant le dernier calcul de P_n(x) = Sigma of y_i * l_i(x)) .

Des suggestions en la rendant plus matlab formelle ?

Grâce

OriginalL'auteur JAN | 2012-06-14

Oui, plusieurs suggestions (mis en œuvre dans la version 1 ci-dessous): if boucle peut être combiné avec for au-dessus d'elle (il suffit de s' index sauter k par quelque chose comme jr(jr~=j) ci-dessous); polynomialSize est toujours égale length(outputConv) qui est toujours égale n (parce que vous avez n points de données, (n-1)th polynôme avec n coefficients), de sorte que le dernier for boucle et la ligne suivante peut également être remplacés par de simples L(k,:) = multiplier * outputConv;

Donc j'ai repris l'exemple sur http://en.wikipedia.org/wiki/Lagrange_polynomial (et adopté leur j-m la notation, mais pour moi j va 1:n et m est 1:n et m~=j), d'où mon initialisation ressemble

clear; clc;
X=[-9 -4 -1 7]; %example taken from http://en.wikipedia.org/wiki/Lagrange_polynomial
Y=[ 5  2 -2 9];
n=length(X); %Lagrange basis polinomials are (n-1)th order, have n coefficients
lj = zeros(1,n); %storage for numerator of Lagrange basis polyns - each w/n coeff
Lj = zeros(n); %matrix of Lagrange basis polyns coeffs (lj(x))
L  = zeros(1,n); %the Lagrange polynomial coefficients (L(x))

alors v 1.0 ressemble

jr=1:n; %j-range: 1<=j<=n
for j=jr %my j is your k
    multiplier = 1;
    outputConv = 1; %numerator of lj(x)
    mr=jr(jr~=j); %m-range: 1<=m<=n, m~=j
    for m = mr %my m is your index
        outputConv = conv(outputConv,[1 -X(m)]);
        multiplier = multiplier * ((X(j) - X(m))^-1);
    end
    Lj(j,:) = multiplier * outputConv; %jth Lagrange basis polinomial lj(x)
end
L = Y*Lj; %coefficients of Lagrange polinomial L(x)

qui peut encore être simplifié si vous vous rendez compte que le numérateur de l_j(x) est un polynôme spécifiques racines - pour qu'il y est une belle commande matlab - poly. De même, le dénominateur est juste que polyn évalué à X(j) - pour qu'il y est polyval. Par conséquent, v 1.9:

jr=1:n; %j-range: 1<=j<=n
for j=jr
    mr=jr(jr~=j); %m-range: 1<=m<=n, m~=j
    lj=poly(X(mr)); %numerator of lj(x)
    mult=1/polyval(lj,X(j)); %denominator of lj(x)
    Lj(j,:) = mult * lj; %jth Lagrange basis polinomial lj(x)
end
L = Y*Lj; %coefficients of Lagrange polinomial L(x)

Pourquoi la version 1.9 et 2.0 pas? eh bien, il y a sans doute une façon de se débarrasser de ce dernier pour la boucle, et d'écrire tout cela en 1 ligne, mais je ne peux pas penser à cela maintenant - c'est une todo pour v 2.0 🙂

Et, pour le dessert, si vous voulez obtenir la même image que wikipédia:

figure(1);clf
x=-10:.1:10;
hold on
plot(x,polyval(Y(1)*Lj(1,:),x),'r','linewidth',2)
plot(x,polyval(Y(2)*Lj(2,:),x),'b','linewidth',2)
plot(x,polyval(Y(3)*Lj(3,:),x),'g','linewidth',2)
plot(x,polyval(Y(4)*Lj(4,:),x),'y','linewidth',2)
plot(x,polyval(L,x),'k','linewidth',2)
plot(X,Y,'ro','linewidth',2,'markersize',10)
hold off
xlim([-10 10])
ylim([-10 10])
set(gca,'XTick',-10:10)
set(gca,'YTick',-10:10)
grid on

produit

La méthode d'interpolation de Lagrange

profiter et n'hésitez pas à réutiliser/améliorer

OriginalL'auteur alexey

Essayer:

X=0:1/20:1; Y=cos(X) and create L and apply polyval(L,1). 
polyval(L,1)=0.917483227909543
cos(1)=0.540302305868140

Pourquoi il y a une énorme différence?

OriginalL'auteur Erdem

Vous devez vous connecter pour publier un commentaire.