LoginSignup
2
1

More than 1 year has passed since last update.

@hibs_MATLAB_Amb(Hibs@MATLAB Ambassador)

ラグランジュ補間で学ぶMATLAB interp

Last updated at Posted at 2022-02-22

はじめに

MATLAB Student Ambassadorとして活動を開始してから1年が経とうとしています. ここまで支えていただいたたくさんの方々に感謝申し上げます. そして, この節目に, これまで私の活動を支えてくれたMATLABコードたちをまとめておこうと考え, 記事にしました.

この記事では, MATLAB R2021bを使用しております.
また, 本記事で示されたコードは, 今回対象とした問題以外での検証等は行っておりません.

ラグランジュ補間

ラグランジュ補間は, 以下の式で定義されています.

f(x) = p_n(x) = \sum_{k=1}^nL_k(x)f_k \\
L_k(x) = \frac{l_k(x)}{l_k(x_k)} \\
l_k(x) = (x-x_1)(x-x_2)\cdots(x-x_{k-1})(x-x_{k+1})\cdots(x-x_n)

ここで, $x_k$は標本点のある位置, そして, $f_k$は$x_k$における標本データを表しています. 今, $x_1\cdots x_k\cdots x_n$という位置において, $f_1\cdots f_k\cdots f_n$という標本データが与えられており, $x'_1\cdots x'_m$における近似値$f'$を計算したいとすると, 上の式で表されている定義は, 下の行列による表現に書き換えることが可能です.

f' = \left[\begin{array}{ccccc}
f_1 & \cdots & f_k & \cdots & f_n
\end{array}\right]
\left[\begin{array}{ccc}
L_1(x'_1) & \cdots & L_1(x'_m) \\
\vdots & \vdots & \vdots \\
L_k(x'_1) & \cdots & L_k(x'_m) \\
\vdots & \vdots & \vdots \\
L_n(x'_1) & \cdots & L_n(x'_m)
\end{array}\right]

これらの式を用いて, MATLABで関数を作成する.

MATLABの補間関数 interp1

MATLABには既にとても使いやすい補間の関数として, interp1が存在します. そこで, 作るラグランジュ補間の関数はできるだけ, interp1に似せた入出力を持たせることにしました (これでMATLABで補間をする方法を学んでくれる人がいると嬉しい).

interp1の基本的な使い方としては,

matlab
vq = interp1(x,v,xq)

のように,

  • $x$ : 標本データのある位置
  • $v$ : 標本データ
  • $xq$ : 近似値の計算が必要な位置

が入力であり,

  • $vq$ : 補間で得られる値

が出力となっています.

作成したMATLAB関数 interpLagVector1

ここでは, $x$, $v$, $xq$がすべてベクトルであると仮定します. この仮定を踏まえて作成したMATLAB関数interpLagVector1が以下になります.

matlab
function vq = interpLagVector1(x,v,xq)

%
% interpLagVector1 - Lagrange Interpolation
%
%   vq = interpLag1(x,v,xq)
%   
%   -Inputs-
%   x   - Sample Points
%   v   - Sample Values
%   xq  - Query Points
%   
%   -Outputs-
%   vq  - Interpolation Value
%
%   --- Usage Example ---
%   
%   x = 0:pi/4:2*pi;
%   v = sin(x);
%   xq = 0:pi/16:2*pi;
%   vq = interpLagVector1(x,v,xq);
%   plot(x,v,'o',xq,vq1,':.');
%
%   ---------------------
%
%

if (isvector(x) && isvector(v) && isvector(xq))
    % Setting
    n = length(x); % Get Number of Cotrol Points
    nq = length(xq); % Get Number of outputting Points
    L = ones(n,nq); % Create Matrix for Lagrange Equation

    % Unify All variables to Row Vector
    x = reshape(x,[1,n]);
    v = reshape(v,[1,n]);
    xq = reshape(xq,[1 nq]);

    % Create Lagrange Equation
    for ii = 1:n
        for jj = 1:nq
            for kk = 1:n
                if (ii ~= kk)
                    L(ii,jj) = L(ii,jj)*(xq(jj)-x(kk))/(x(ii)-x(kk));
                end
            end
        end
    end

    vq = v*L; % Get Answer
else
    error('Error!! All inputs should be Vector.');
end
end

ここで, MATLABでは, 関数の作成と同時に, その中身の説明を加えることができます.

interpLagVector1の確認

ここでは, ベッセル関数(besselj)を用いて作成したinterpLagVector1を試してみました.

matlab
%% Setting
start = 1;
stop = 5;

%% Single Data Interpolation 
% Create Sampling Data
x = start:stop;
v = besselj(1,x);

% Exact Solution
x_ex = start:0.1:stop;
v_ex = besselj(1,x_ex);

% Interpolation
xq = start:0.25:stop;
vq = interpLagVector1(x,v,xq);

% Visualize
figure
hold on
plot(x_ex,v_ex,'k-','LineWidth',1.3);
plot(x,v,'bo','MarkerSize',10);
plot(xq,vq,'rx','LineWidth',2.0,'MarkerSize',7);
hold off
grid on
title('Lagrange Interpolation');
xlabel('x');
ylabel('f');
legend("Exact Solution","Sampling Data","Interpolation Data");

qiita_lagrange_example_single.png

MATLAB補間関数への行列の入力

MATLABにある補間の関数interp1では, $v$及び$xq$について行列を入力することも可能です. この場合, 複数のデータの補間を計算したり, あるデータから異なる点での補間の値を算出することができます. そこで, 今回は, MATrix LABoratoryで行列を入力としたラグランジュ補間を行ってみました.

作成した関数 interpLag1

ここでは, $x$はベクトルであり, $v$, $xq$のうち少なくとも片方はベクトルであると仮定します. この仮定を踏まえて作成したMATLAB関数interpLag1が以下になります. この関数では先に作成したinterpLagVector1を用いています.

matlab
function vq = interpLag1(x,v,xq)

%
% interpLag1 - Lagrange Interpolation
%
%   vq = interpLag1(x,v,xq)
%   
%   -Inputs-
%   x   - Sample Points
%   v   - Sample Values
%   xq  - Query Points
%   
%   -Outputs-
%   vq  - Interpolation Value
%
%   --- Usage Example ---
%   
%   x = 0:pi/4:2*pi;
%   v = sin(x);
%   xq = 0:pi/16:2*pi;
%   vq = interpLag1(x,v,xq);
%   plot(x,v,'o',xq,vq1,':.');
%
%   ---------------------
%
%

if (isvector(x))
    if (length(x) == length(unique(x)))
        if (isvector(v) && isvector(xq))
            vq = interpLagVector1(x,real(v),xq) + 1i*interpLagVector1(x,imag(v),xq);
        elseif (isvector(v))
            vq = zeros(size(xq));
            for ll=1:length(xq(1,:))
                vq(:,ll) = interpLagVector1(x,real(v),xq(:,ll)) + 1i*interpLagVector1(x,imag(v),xq(:,ll));
            end
        elseif (isvector(xq))
            vq = zeros(length(xq),length(v(1,:)));
            for ll=1:length(v(1,:))
                vq(:,ll) = interpLagVector1(x,real(v(:,ll)),xq) + 1i*interpLagVector1(x,imag(v(:,ll)),xq);
            end
        else
            error('Error!! Either v or xq should be vector.');
        end

    else
        error('Error!! x should not has same data.');
    end
else
    error('Error!! The input should be Vector.');
end
end

% Lagrange Interpolation
function vq = interpLagVector1(x,v,xq)
if (isvector(x) && isvector(v) && isvector(xq))
    % Setting
    n = length(x); % Get Number of Cotrol Points
    nq = length(xq); % Get Number of outputting Points
    L = ones(n,nq); % Create Matrix for Lagrange Equation

    % Unify All variables to Row Vector
    x = reshape(x,[1,n]);
    v = reshape(v,[1,n]);
    xq = reshape(xq,[1 nq]);

    % Create Lagrange Equation
    for ii = 1:n
        for jj = 1:nq
            for kk = 1:n
                if (ii ~= kk)
                    L(ii,jj) = L(ii,jj)*(xq(jj)-x(kk))/(x(ii)-x(kk));
                end
            end
        end
    end

    vq = v*L; % Get Answer
else
    error('Error!! All inputs should be Vector.');
end
end

interpLag1の確認

ここでも, ベッセル関数(besselj)を用いて, 作成したinterpLag1を試してみました.

複数のデータのラグランジュ補間

ここでは, 1次, 2次及び3次のベッセル関数から標本データ$v$を取得し, ベクトル$xq$の各点における補間の値$vq$を算出しました.

matlab
clear all; close all; clc;
%% Setting
start = 1;
stop = 5;

%% Multiple Sample Interpolation
% Number of Data
n = 3;

% Create Sampling Data
x = start:stop;
v = zeros(length(x),n);
for ii = 1:n
    v(:,ii) = besselj(ii,x);
end

% Exact Solution
x_ex = start:0.1:stop;
v_ex = zeros(length(x_ex),n);
for ii=1:n
    v_ex(:,ii) = besselj(ii,x_ex);
end

% Interpolation
xq = start:0.25:stop;
vq = interpLag1(x,v,xq);

% Visualize
figure
hold on
plot(x_ex,v_ex,'-','LineWidth',1.7);
plot(x,v,'o','MarkerSize',10);
plot(xq,vq,'x','LineWidth',2.0,'MarkerSize',7);
hold off
grid on
title('Lagrange Interpolation');
xlabel('x');
ylabel('f');

qiita_lagrange_example_multipleSample.png

図において, 丸が標本データ$v$, xが補間で得られた値$vq(:,1), vq(:,2), vq(:,3)$, 実践が1次, 2次及び3次のベッセル関数を示しています.

あるデータから複数の点の組み合わせでの補間

ここでは, 先ほどまでと同様に, $1 \leq x \leq 5$の範囲を等分割した列ベクトル$xq(:,1)$とランダムに同数の点を抽出した列ベクトル$xq(:,2)$を作成しました. そして, この2種類のベクトルの各点における補間した値$vq$を算出しました.

matlab
clear all; close all; clc;
%% Setting
start = 1;
stop = 5;

%% Get Multiple Interpolation Data from Single Sample
% Create Sampling Data
x = start:stop;
v = besselj(1,x);

% Exact Solution
x_ex = start:0.1:stop;
v_ex = besselj(1,x_ex);

% Interpolation
xq1 = start:0.25:stop;
xq2 = sort(start+(stop-start)*rand(1,length(xq1)));
xq = [xq1' xq2'];
vq = interpLag1(x,v,xq);

% Visualize
figure
subplot(2,1,1)
hold on
plot(x_ex,v_ex,'k-','LineWidth',1.3);
plot(x,v,'bo','MarkerSize',10);
plot(xq(:,1),vq(:,1),'rx','LineWidth',2.0,'MarkerSize',7);
hold off
grid on
title('Lagrange Interpolation');
xlabel('x');
ylabel('f');
legend("Exact Solution","Sampling Data","Interpolation Data 1");

subplot(2,1,2)
hold on
plot(x_ex,v_ex,'k-','LineWidth',1.3);
plot(x,v,'bo','MarkerSize',10);
plot(xq(:,2),vq(:,2),'rx','LineWidth',2.0,'MarkerSize',7);
hold off
grid on
title('Lagrange Interpolation');
xlabel('x');
ylabel('f');
legend("Exact Solution","Sampling Data","Interpolation Data 2");

qiita_lagrange_example_multipleInt.png

図において, 上側は$1 \leq x \leq 5$を等分割した点での補間により得られた値$vq(:,1)$, 下側は$1 \leq x \leq 5$からランダムに抽出した点での補間により得られた値$vq(:,2)$を示している.

おわりに

今回は, 「$v$または$xq$の少なくとも一方がベクトルである」という仮定をつけてしまいましたが, 本家interp1ではn次元配列も使用することが可能です (難しかったので妥協してしまいました).

私を支えてくれている関数はまだまだたくさんあるので, 少しずつ備忘録として残しておこうと思います.
また, 今回は, 「MATLABにラグランジュ補間の関数がないなー」と思って作成しましたが, もしあれば教えていただきたいです!!

最後に, 本記事のアイデアをくださり, 一緒に考えてくれた方に, 心より感謝申し上げます.

参考文献

E. クライツィグ著, 田村義保訳, 技術者のための高等数学=5 数値解析 原書第8版, 培風館, 2015
MATLAB Documentation interp1

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
What you can do with signing up
Sign upLogin
2
1