単調な三次補間

数値解析の数学的分野において、単調 3 次補間は、補間される データ セットの単調性を維持する3 次補間の変形です。

単調性は線形補間によって維持されますが、三次補間では保証されません。

単調補間は、結果として得られるエルミート スプラインの単調性を保証するために 接線を修正した3 次エルミート スプラインを使用して実現できます。${\displaystyle m_{i}}$

単調五次エルミート補間のためのアルゴリズムも利用できます。

各データ ポイントの補間接線を選択する方法はいくつかあります。このセクションでは、Fritsch-Carlson 法の使用法について概説します。アルゴリズムは 1 回だけ実行する必要があることに注意してください。

データ ポイントをソート順にインデックス付けします。 ${\displaystyle (x_{k},y_{k})}$${\displaystyle k=1,\,\dots \,n}$

1. 連続する点間の割線の傾きを計算します。

   ${\displaystyle \delta_{k}={\frac {y_{k+1}-y_{k}}{x_{k+1}-x_{k}}}}$

   のために。${\displaystyle k=1,\,\dots \,n-1}$
   
   
2. これらの割り当ては暫定的なものであり、残りの手順で置き換えられる可能性があります。すべての内部データポイントの接線をセカントの平均として初期化します。

   ${\displaystyle m_{k}={\frac {\delta _{k-1}+\delta _{k}}{2}}}$

   エンドポイントについては、片側差分を使用します。${\displaystyle k=2,\,\dots \,n-1}$
   
   

   ${\displaystyle m_{1}=\delta _{1}\quad {\text{ および }}\quad m_{n}=\delta _{n-1}\,}$。

   と が異符号の場合はと設定します。${\displaystyle \delta_{k-1}}$${\displaystyle \delta_{k}}$${\displaystyle m_{k}=0}$
   
   
3. については、 （連続する 2 つの点が等しい場合）常に、これらの点を結ぶスプラインが単調性を保つために平坦でなければならないため、 と設定されます。これらの については、手順 4 と 5 は無視してください。${\displaystyle k=1,\,\dots \,n-1}$${\displaystyle \delta_{k}=0}$${\displaystyle y_{k}=y_{k+1}}$
   ${\displaystyle m_{k}=m_{k+1}=0,}$
   ${\displaystyle k\,}$
   
   
4. させて

   ${\displaystyle \alpha _{k}=m_{k}/\delta _{k}\quad {\text{ および }}\quad \beta _{k}=m_{k+1}/\delta _{k}}$。

   または が負の場合、入力データ ポイントは厳密には単調ではなく、 は局所的極値になります。このような場合でも、の場合またはの場合を選択することで、区分的単調曲線を生成できますが、厳密な単調性は全体的には不可能です。${\displaystyle \alpha_{k}}$${\displaystyle \beta _{k}}$${\displaystyle (x_{k},\,y_{k})}$${\displaystyle m_{k}=0\,}$${\displaystyle \alpha _{k}$${\displaystyle m_{k+1}=0\,}$${\displaystyle \beta _{k}$
   
   
5. オーバーシュートを防ぎ、単調性を確保するには、次の 3 つの条件のうち少なくとも 1 つを満たす必要があります。

（a）関数

${\displaystyle \phi _{k}=\alpha _{k}-{\frac {(2\alpha _{k}+\beta _{k}-3)^{2}}{3(\alpha _ {k}+\beta _{k}-2)}}>0\,}$、 または

（ b）または ${\displaystyle \alpha _{k}+2\beta _{k}-3\leq 0\,}$

（ハ）.${\displaystyle 2\alpha _{k}+\beta _{k}-3\leq 0\,}$

厳密な単調性を保証するには、条件 (a) のみが正である必要があります。${\displaystyle \phi_{k}}$

この制約を満たす簡単な方法は、ベクトルを半径3の円に制限することです。つまり、 の場合、次のように設定します。${\displaystyle (\alpha _{k},\,\beta _{k})}$${\displaystyle \alpha _{k}^{2}+\beta _{k}^{2}>9\,}$

${\displaystyle \tau _{k}={\frac {3}{\sqrt {\alpha _{k}^{2}+\beta _{k}^{2}}}}\,}$、

そして接線を次のように再スケールする。

${\displaystyle m_{k}=\tau _{k}\,\alpha _{k}\,\delta _{k}\quad {\text{ and }}\quad m_{k+1}=\tau _{k}\,\beta _{k}\,\delta _{k}\,}$。

あるいは、と を制限すれば十分です。これを実現するには、の場合は を設定し、の場合は を設定します。${\displaystyle \alpha _{k}\leq 3}$${\displaystyle \beta _{k}\leq 3\,}$${\displaystyle \alpha _{k}>3\,}$${\displaystyle m_{k}=3\,\delta _{k}\,}$${\displaystyle \beta _{k}>3\,}$${\displaystyle m_{k+1}=3\,\delta _{k}\,}$

上記の前処理の後、補間されたスプラインの評価は、、、 およびのデータを使用して、 3次エルミートスプラインと同等になります。 ${\displaystyle x_{k}}$${\displaystyle y_{k}}$${\displaystyle m_{k}}$${\displaystyle k=1,\,\dots \,n}$

で を評価するには、 、 、 の間にあるシーケンス内のインデックス、つまり を見つけます。計算します。 ${\displaystyle x}$${\displaystyle k}$${\displaystyle x}$${\displaystyle x_{k}}$${\displaystyle x_{k+1}}$${\displaystyle x_{k}\leq x\leq x_{k+1}}$

${\displaystyle \Delta =x_{k+1}-x_{k}\quad {\text{ and }}\quad t={\frac {x-x_{k}}{\Delta }}}$

補間値は

${\displaystyle f_{\text{interpolated}}(x)=y_{k}\cdot h_{00}(t)+\Delta \cdot m_{k}\cdot h_{10}(t)+y_{k+1}\cdot h_{01}(t)+\Delta \cdot m_{k+1}\cdot h_{11}(t)}$

ここで、3次エルミートスプラインの基底関数です。 ${\displaystyle h_{ii}}$

次のJavaScript実装は、データ セットを受け取り、単調な 3 次スプライン補間関数を生成します。

/* 
* 単調 3 次スプライン補間
* 使用例は下部に記載されています。これは完全に機能するパッケージです。 * たとえば、これは
https://www.programiz.com/javascript/online-compiler/などのサイトで実行するか、nodeJS を使用して実行できます。 
*/ function DEBUG ( s ) { /* ソルバーの詳細な出力を有効にするには、次のコメントを解除します: */ //console.log(s); } var j = 0 ; var createInterpolant = function ( xs , ys ) { var i , length = xs . length ; // 長さの問題に対処しますif ( length != ys . length ) { throw 'xs と ys の数が等しくなければなりません。' ; } if ( length === 0 ) { return function ( x ) { return 0 ; }; } if ( length === 1 ) { // 実装: 結果を事前計算することで、ys が後で変更された場合の問題を防ぎ、ys のガベージコレクションを可能にします// 実装: 単項加算により、値を適切に数値に変換しますvar result = + ys [ 0 ]; return function ( x ) { return result ; }; } // xs と ys を並べ替えて、xs がソートされるようにしますvar indexes = []; for ( i = 0 ; i < length ; i ++ ) { indexes . push ( i ); } indexes . sort ( function ( a , b ) { return xs [ a ] < xs [ b ] ? - 1 : 1 ; }); var oldXs = xs


  
    
    

   
     
        
    
    
           
               
        
        
        
           
             
    
    
    
       
              
               
       , oldYs = ys ; // 実装: 新しい配列を作成すると、入力配列が後で変更された場合の問題も防止されますxs = []; ys = []; // 実装: 単項加算により、値が数値に適切に変換されますfor ( i = 0 ; i < length ; i ++ ) { xs . push ( + oldXs [ indexes [ i ]]); ys . push ( + oldYs [ indexes [ i ]]); }   
    
         
    
               

    DEBUG ( "debug: xs = [ " + xs + " ]" ) DEBUG ( "debug: ys = [ " + ys + " ]" ) // 連続した差と傾きを取得しますvar dys = [], dxs = [], ms = []; for ( i = 0 ; i < length - 1 ; i ++ ) { var dx = xs [ i + 1 ] - xs [ i ], dy = ys [ i + 1 ] - ys [ i ]; dxs . push ( dx ); dys . push ( dy ); ms . push ( dy / dx ); } // 次数 1 の係数を取得しますvar c1s = [ ms [ 0 ]]; for ( i = 0 ; i < dxs . length - 1 ; i ++ ) { var m = ms [ i ], mNext = ms [ i + 1 ]; if ( m * mNext <= 0 ) { c1s . push ( 0 ); } else { var dx_ = dxs [ i ], dxNext = dxs [ i + 1 ], common = dx_ + dxNext ; c1s . push ( 3 * common / (( common + dxNext ) / m + ( common + dx_ ) / mNext    
        
    
    
             
              
                      
          
    
    
       
              
                
            
            
          
                         
                  
) ) ; } } c1s.push ( ms [ ms.length - 1 ] ) ;        
    
      

    DEBUG ( "debug: dxs = [ " + dxs + " ]" ) DEBUG ( "debug: ms = [ " + ms + " ]" ) DEBUG ( "debug: c1s.length = " + c1s . length ) DEBUG ( "debug: c1s = [ " + c1s + " ]" ) // 2 次係数と 3 次係数を取得しますvar c2s = [], c3s = []; for ( i = 0 ; i < c1s . length - 1 ; i ++ ) { var c1 = c1s [ i ]; var m_ = ms [ i ]; var invDx = 1 / dxs [ i ]; var common_ = c1 + c1s [ i + 1 ] - m_ - m_ ; DEBUG ( "debug: " + i + ". c1 = " + c1 ); DEBUG ( "debug: " + i + ". m_ = " + m_ ); DEBUG ( "debug: " + i + ". invDx = " + invDx ); DEBUG ( "debug: " + i + ". common_ = " + common_ ); c2s . push (( m_ - c1 - common_ ) * invDx ); c3s . push ( common_ * invDx * invDx ); } DEBUG ( "debug: c2s = [ " + c2s + " ]" ) DEBUG ( "debug: c3s = [ " + c3s + " ]" )    
        
      
        
    
    
          
              
           
           
           
                   
              
              
              
              
            
        
    
        
        

    // 補間関数を返します
return function ( x ) { // データ セットの右端のポイントは正確な結果を返しますvar i = xs . length - 1 ; //if (x == xs[i]) { return ys[i]; } // x が含まれる区間を検索し、x が元の xs の 1 つである場合は対応する y を返しますvar low = 0 , mid , high = c3s . length - 1 , rval , dval ; while ( low <= high ) { mid = Math . floor ( 0.5 * ( low + high )); var xHere = xs [ mid ]; if ( xHere < x ) { low = mid + 1 ; } else if ( xHere > x ) { high = mid - 1 ; } else { j ++ ; i = mid ; var diff = x - xs [ i ]; rval = ys [ i ] + diff * ( c1s [ i ] + diff * ( c2s [ i ] + diff * c3s [ i ])); dval = c1s [ i ] + diff * ( 2 * c2s [ i ] + diff * 3 * c3s [ i ]); DEBUG ( "debug: " + j + ". x = " + x + ". i = " + i + ", diff = " + diff +      
        
             
        
        
        
                   
            
                
               
                      
                       
             
                
                  
                     
                               
                          
                                "、rval = " + rval + "、dval = " + dval ); return [ rval , dval ]; i =数学。​​最大( 0 、高);      
                    
            
        
           

        // 補間
var diff = x - xs [ i ]; j ++ ; rval = ys [ i ] + diff * ( c1s [ i ] + diff * ( c2s [ i ] + diff * c3s [ i ])); dval = c1s [ i ] + diff * ( 2 * c2s [ i ] + diff * 3 * c3s [ i ]); DEBUG ( "debug: " + j + ". x = " + x + ". i = " + i + ", diff = " + diff + ", rval = " + rval + ", dval = " + dval ); return [ rval , dval ]; }; };             
        
                       
                  
                              
            
    


/*
   以下の使用例では、0 <= x <= 4 の場合に x^2 を近似します。

   コマンドラインの使用例（Node.js のインストールが必要です）: 
   node monotone-cubic-spline.js 
*/

var X = [ 0 , 1 , 2 , 3 , 4 ]; var F = [ 0 , 1 , 4 , 9 , 16 ] ; var f = createInterpolant ( X , F ); var N = X.length ; console.log ( "# BLOCK 0 :: monotone-cubic-spline.js のデータ"); console.log ( " X " + " \ t " + "F" ); for ( var i = 0 ; i < N ; i + = 1 ) { console.log ( F [ i ] + ' \t' + X[i]); } console.log ( " " ) ; console.log ( " " ) ; console.log ( " # BLOCK 1 :: monotone - cubic- spline.jsの補間データ" ) ; console .log ( " x " + "\t\t" + " P(x) " + "\t\t" + " dP(x)/dx " ); var message = '' ; var M = 25 ; for ( var i = 0 ; i <= M ; i += 1 ) { var x = X [ 0 ] + ( X [ N - 1 ] - X [ 0 ]) * i / M ; var rvals = f ( x ); var P = rvals [ 0 ]; var D       
       
   
   

    
           
        




        
   
   
           
         
       
       
      = rvals [ 1 ]; message += x .toPrecision ( 15 ) + '\t' + P .toPrecision ( 15 ) + ' \ t ' + D .toPrecision ( 15 ) + ' \n' ; } console.log ( message ) ; 
                

• Fritsch, FN; Carlson, RE (1980). 「単調区分三次補間」. SIAM Journal on Numerical Analysis . 17 (2). SIAM: 238–246. doi :10.1137/0717021.
• Dougherty , RL; Edelman, A.; Hyman, JM (1989 年 4 月)。「正値性、単調性、または凸性を維持する 3 次および 5 次エルミート補間」。計算数学。52 (186): 471–494。doi : 10.2307/2008477。

• GPLv 2 ライセンスのC++実装: MonotCubicInterpolator.cpp MonotCubicInterpolator.hpp