はじめに

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");

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');

図において, 丸が標本データ$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");

図において, 上側は$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

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

はじめに

ラグランジュ補間

MATLABの補間関数 interp1

作成したMATLAB関数 interpLagVector1

interpLagVector1の確認

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

作成した関数 interpLag1

interpLag1の確認

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

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

おわりに

参考文献