エクセルマクロ(Excel VBA)でラグランジュ補間を行う

#エクセルマクロ(Excel VBA)でラグランジュ補間を行う

はじめに

エクセルでラグランジュ補間を行う必要があったので、マクロを作成しました。
関数を指定し、求めたい地点とデータを指定するだけで簡単に補間が行えます。
ラグランジュ補間は点数を多くすると解の振動が発生するので、標準では近傍4点のみで補間を行うようにしています。

関数の仕様

関数名：LagrangeInterpolation(x, 既知のy, 既知のx, 近傍のみで補間するか)
引数：
x: 求めたい地点
既知のy: データの範囲(Y)
既知のx: データの範囲(X)
近傍のみ: Falseにすると全ての地点で補間する。指定しないと自動的にTrue

※引数が"x""既知のy""既知のx" の順番になっているのはFORECAST関数と同一の引数にしたいからです。

ソース

Lagrange補間.bas

Option Explicit

Function LagrangeInterpolation(X As Double, 既知のy As Range, 既知のx As Range, Optional 近傍のみ As Boolean = True) As Variant
    
    ' 近傍のみで補間する場合は前後N/2点ずつ補間する
    ' 例：N=4の場合は前後2点ずつ
    ' Nは2の倍数である必要がある
    Const N = 4
      
    ' yとxのデータは同じ個数であることが前提
    If 既知のx.count <> 既知のy.count Then
        ' xとyでデータ個数が異なる場合はN/Aを返す
        LagrangeInterpolation = CVErr(xlErrNA)
        Exit Function
    End If
    
    ' 生データの個数を代入
    Dim DataNum As Long
    DataNum = 既知のy.count
    
    ' ループ関数を定義
    Dim i As Long
    Dim j As Long
    
    ' Rangeデータを配列に保存
    Dim Xdatas() As Double
    ReDim Xdatas(既知のx.count)
    Dim Ydatas() As Double
    ReDim Ydatas(既知のy.count)
    For i = 1 To 既知のx.count
        Xdatas(i) = 既知のx(i)
        Ydatas(i) = 既知のy(i)
    Next
    
    ' 地点をXでソートし、X前後の地点をN点絞り込む
    ' 地点数Nよりデータ数のほうが少ない場合は何もしない
    If 近傍のみ And N <= DataNum Then
        
        ' ソートを行う
        Dim BufDouble As Double
        For i = 1 To UBound(Xdatas)
            For j = 1 To UBound(Xdatas)
                If Xdatas(i) < Xdatas(j) Then
                    BufDouble = Xdatas(i)
                    Xdatas(i) = Xdatas(j)
                    Xdatas(j) = BufDouble
                    
                    BufDouble = Ydatas(i)
                    Ydatas(i) = Ydatas(j)
                    Ydatas(j) = BufDouble
                End If
            Next
        Next
        
        ' 切り出した後のXとYのデータ
        Dim CuttedXDatas() As Double
        ReDim CuttedXDatas(N)
        Dim CuttedYDatas() As Double
        ReDim CuttedYDatas(N)
        
        ' Xに最も近い要素番号を求める
        Dim StartCutNum As Long
        If X < Xdatas(2) Then
            ' Xより小さいデータが十分に存在しない場合、1～Nまでで補間を行う
            StartCutNum = 1
        ElseIf Xdatas(UBound(Xdatas) - 2) < X Then
            ' Xより大きいデータが十分に存在しない場合、生データ数-N～生データ数までで補間を行う
            StartCutNum = UBound(Xdatas) - N + 1
        Else
            ' X近傍に十分なデータ数がある場合、前後で補間
            For i = 1 To UBound(Xdatas) - 1
                If Xdatas(i) < X And X <= Xdatas(i + 1) Then
                    StartCutNum = i
                    Exit For
                End If
            Next
        End If
        
        If StartCutNum = 1 Then
            ' 1～N点で補間を行う
            For i = 1 To N
                CuttedXDatas(i) = Xdatas(i)
                CuttedYDatas(i) = Ydatas(i)
            Next
        ElseIf StartCutNum = UBound(Xdatas) - N + 1 Then
            ' 生データ数-N～生データ数で補間を行う
            For i = StartCutNum To DataNum
                CuttedXDatas(i - StartCutNum + 1) = Xdatas(i)
                CuttedYDatas(i - StartCutNum + 1) = Ydatas(i)
            Next
        Else
            ' 前後の地点で補間を行う
            For i = StartCutNum - N / 2 To StartCutNum + 1
                CuttedXDatas(i - StartCutNum + N / 2 + 1) = Xdatas(i)
                CuttedYDatas(i - StartCutNum + N / 2 + 1) = Ydatas(i)
            Next
        End If
        
        ' 切り出したデータを反映
        DataNum = N
        ReDim Xdatas(N)
        Xdatas = CuttedXDatas
        ReDim Ydatas(N)
        Ydatas = CuttedYDatas
        
    End If
    
    ' ラグランジュ既定関数の分子および分母
    Dim Numerator As Double     ' 分子
    Dim Denominator As Double   ' 分母
    
    ' 関数の戻り値を初期化
    LagrangeInterpolation = 0
    
    ' ラグランジュ既定関数を計算し、yの値との積を合計する
    For i = 1 To DataNum
        Numerator = 1
        Denominator = 1
        For j = 1 To DataNum
            If i <> j Then
                Numerator = Numerator * (X - Xdatas(j))
                Denominator = Denominator * (Xdatas(i) - Xdatas(j))
            End If
        Next
        LagrangeInterpolation = LagrangeInterpolation + (Ydatas(i) * Numerator / Denominator)
    Next
    
End Function

使用例および結果

この関数を用いて、以下の数式を補間してみましょう。

\frac{1}{1+8x^2}

###エクセルのスクリーンショット

補間後のグラフ

黒点線：真値
薄赤太線：すべての点で補間した結果（解が振動する）
薄青太線：近傍４点で補間した結果（解が振動しない）

まとめ

データ数が多いと、ただ単にラグランジュ補間を行うだけでは解が振動してしまいます。
そのため今回は標準で、近傍地点のみで補間をするように工夫しました。
ただし、近傍地点のみで補間しても振動を避けられないこともあるので、解の取り扱いには注意が必要です。

2023/8/20追記

セル内にNAなどのエラーが含まれているときにその値を無視するバージョンを作成したので共有します。

ソース

Lagrange補間.bas

Option Explicit

Function LagrangeInterpolation(X As Double, 既知のy As Range, 既知のx As Range, Optional 近傍のみ As Boolean = True) As Variant

    ' 近傍のみで補間する場合は前後N/2点ずつ補間する
    ' 例：N=4の場合は前後2点ずつ
    ' Nは2の倍数である必要がある
    Const N = 4

    ' yとxのデータは同じ個数であることが前提
    If 既知のx.Count <> 既知のy.Count Then
        ' xとyでデータ個数が異なる場合はN/Aを返す
        LagrangeInterpolation = CVErr(xlErrNA)
        Exit Function
    End If

    ' 生データの個数を代入
    Dim DataNum As Long
    DataNum = 既知のy.Count

    ' ループ関数を定義
    Dim i As Long
    Dim j As Long

    ' Rangeデータを配列に保存
    Dim Xdatas() As Double
    ReDim Xdatas(既知のx.Count)
    Dim Ydatas() As Double
    ReDim Ydatas(既知のy.Count)
    For i = 1 To 既知のx.Count
        If Not (IsNumeric(既知のx(i)) And IsNumeric(既知のy(i))) Then
            GoTo Continue
        End If
        
        Xdatas(i) = 既知のx(i)
        Ydatas(i) = 既知のy(i)
Continue:
    Next

    ' 地点をXでソートし、X前後の地点をN点絞り込む
    ' 地点数Nよりデータ数のほうが少ない場合は何もしない
    If 近傍のみ And N <= DataNum Then

        ' ソートを行う
        Dim BufDouble As Double
        For i = 1 To UBound(Xdatas)
            For j = 1 To UBound(Xdatas)
                If Xdatas(i) < Xdatas(j) Then
                    BufDouble = Xdatas(i)
                    Xdatas(i) = Xdatas(j)
                    Xdatas(j) = BufDouble

                    BufDouble = Ydatas(i)
                    Ydatas(i) = Ydatas(j)
                    Ydatas(j) = BufDouble
                End If
            Next
        Next

        ' 切り出した後のXとYのデータ
        Dim CuttedXDatas() As Double
        ReDim CuttedXDatas(N)
        Dim CuttedYDatas() As Double
        ReDim CuttedYDatas(N)

        ' Xに最も近い要素番号を求める
        Dim StartCutNum As Long
        If X < Xdatas(2) Then
            ' Xより小さいデータが十分に存在しない場合、1～Nまでで補間を行う
            StartCutNum = 1
        ElseIf Xdatas(UBound(Xdatas) - 2) < X Then
            ' Xより大きいデータが十分に存在しない場合、生データ数-N～生データ数までで補間を行う
            StartCutNum = UBound(Xdatas) - N + 1
        Else
            ' X近傍に十分なデータ数がある場合、前後で補間
            For i = 1 To UBound(Xdatas) - 1
                If Xdatas(i) < X And X <= Xdatas(i + 1) Then
                    StartCutNum = i
                    Exit For
                End If
            Next
        End If

        If StartCutNum = 1 Then
            ' 1～N点で補間を行う
            For i = 1 To N
                CuttedXDatas(i) = Xdatas(i)
                CuttedYDatas(i) = Ydatas(i)
            Next
        ElseIf StartCutNum = UBound(Xdatas) - N + 1 Then
            ' 生データ数-N～生データ数で補間を行う
            For i = StartCutNum To DataNum
                CuttedXDatas(i - StartCutNum + 1) = Xdatas(i)
                CuttedYDatas(i - StartCutNum + 1) = Ydatas(i)
            Next
        Else
            ' 前後の地点で補間を行う
            For i = StartCutNum - N / 2 To StartCutNum + 1
                CuttedXDatas(i - StartCutNum + N / 2 + 1) = Xdatas(i)
                CuttedYDatas(i - StartCutNum + N / 2 + 1) = Ydatas(i)
            Next
        End If

        ' 切り出したデータを反映
        DataNum = N
        ReDim Xdatas(N)
        Xdatas = CuttedXDatas
        ReDim Ydatas(N)
        Ydatas = CuttedYDatas

    End If

    ' ラグランジュ既定関数の分子および分母
    Dim Numerator As Double     ' 分子
    Dim Denominator As Double   ' 分母

    ' 関数の戻り値を初期化
    LagrangeInterpolation = 0

    ' ラグランジュ既定関数を計算し、yの値との積を合計する
    For i = 1 To DataNum
        Numerator = 1
        Denominator = 1
        For j = 1 To DataNum
            If i <> j Then
                Numerator = Numerator * (X - Xdatas(j))
                Denominator = Denominator * (Xdatas(i) - Xdatas(j))
            End If
        Next
        LagrangeInterpolation = LagrangeInterpolation + (Ydatas(i) * Numerator / Denominator)
    Next

End Function