% Lagrange 補間の動作サンプル
clear all;
clc;
close all;


% 諸パラメータ（ここをいろいろ変更してみるとよい）
n = 10;         % 標本点数；n-1次多項式が得られる

% 標本点の設定：どちらかを選択してみよう
xx = linspace(-1, 1, n);    % 等間隔点
%xx = cos( linspace(1, 2*n-1, n) * pi / (2 * n) );   % Chebyshev 点

% データの生成：被近似関数を選んで様子を見よう
%func = @func1;     % 易しい
func = @func2;     % 難しい
yy = func(xx);      % データ生成

% Lagrange 補間多項式の生成
%   【参考】MATLAB における多項式表現の持ち方
%           a_0 x^n + a_1 x^{n-1} + ... + a_n x^0 = 0 を [a_0 a_1 ... a_n]
%           というベクトルで持つ．
%           poly(ベクトル)：「ベクトル」要素を根に持つ多項式（a_0 = 1) を生成
%           poyyval：上の表現形式の多項式を指定の点で評価する
p = 0;          % MATLAB の多項式表現を入れる器
for i = 1:n
    xl = xx; xl(i) = [];    % 標本点リストから i 番目の点だけ抜く
    L = poly(xl);           % そのリストを根に持つ
    for j = [1:i-1, i+1:n]  % k を除いて for ループを回す（こういう書き方ができる）
        L = L / (xx(i) - xx(j));    % Lagrange多項式に係数整形
                            % MATLAB はループを回したら負けなのでもっと言いやり方を考えたいところ
    end
    p = p + yy (i) * L;     % l_i(x) を追加
end

% 補間多項式と厳密会の図示
xxx = linspace(-1, 1, 1000);     % 図示に使う点（滑らかに見せるためたくさん取る）
yt = func(xxx);                  % 被近似関数の値
yp = polyval(p, xxx);            % 補間多項式の値（Horner法で評価）
plot(xxx, yt,'b');               % 被近似関数をプロット
hold on;
plot(xxx, yp,'r');               % Lagrange補間のプロット
plot(xx, yy, 'bo');              % 補間点を表示
legend('被近似関数','Lagrange補間','補間点');
title(['標本点数 = ', num2str(n)]);
xlabel('x');
ylabel('y');
ylim([-0.4 1.2]);                % y 方向の表示範囲を限定（＝図のスケールを揃える）


% 以降，途中で使われるサブルーチン：被近似関数の定義
function y = func1(x)
    y = 1.0 ./ (1.0 + x.^2);    % y = 1/(1+x^2) （易しい）
end

function y = func2(x)
    y = 1.0 ./ (1.0 + 25.0 * x.^2);    % y = 1/(1+25x^2)　（難しい）
end