Monotone kubisk interpolation

Inom det matematiska området för numerisk analys är monoton kubisk interpolation en variant av kubisk interpolation som bevarar monotoniteten hos datamängden som interpoleras .

Monotonicitet bevaras genom linjär interpolation men garanteras inte av kubisk interpolation .

Monotone kubisk Hermite interpolation

Exempel som visar icke-monoton kubisk interpolation (i rött) och monoton kubisk interpolation (i blått) av en monoton datamängd.

Monotone interpolation kan åstadkommas med kubisk Hermite spline med tangenterna $m_{i}$ modifierade för att säkerställa monotoniteten hos den resulterande Hermite spline.

En algoritm är också tillgänglig för monoton kvintisk Hermite-interpolation.

Interpolant val

Det finns flera sätt att välja interpolerande tangenter för varje datapunkt. Detta avsnitt kommer att beskriva användningen av Fritsch–Carlson-metoden. Observera att endast ett pass av algoritmen krävs.

Låt datapunkterna vara $(x_{k},y_{k})$ indexerade i sorterad ordning för ${\displaystyle k=1,\,\dots$ .

1. Beräkna lutningarna för sekantlinjerna mellan på varandra följande punkter:

   $\delta _{k}={\frac {y_{k+1}-y_ {k}}{x_{k+1}-x_{k}}}$

   för $k=1,\,\dots \,n-1$ .
   
   
2. 
   
   Dessa tilldelningar är provisoriska och kan ersättas i de återstående stegen. Initiera tangenterna vid varje inre datapunkt som medelvärdet av sekanterna,

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

   för $k=2,\,\dots \,n-1$ . För endpoints, använd ensidiga skillnader:

   $m_{1}=\delta _{1}\quad {\text{ och }}\quad m_{ n}=\delta _{n-1}\,$ .

   Om $\delta _{k-1}$ och $\delta _{k}$ har motsatta tecken, sätt $m_{k}=0$ .
   
   
3. 
   
   För $k=1,\,\dots \,n-1$ , där någonsin $\delta _{k}=0$ (där alltid två på varandra följande $y_{k}=y_{k+1}$ är lika), set $m_{k}=m_{k+1}= 0,$ eftersom spline som förbinder dessa punkter måste vara platt för att bevara monotoniteten. Ignorera steg 4 och 5 för de $k\,$ .
   
   
4. Låt

   $\alpha _{k}=m_{k}/\delta _{k}\quad {\text{ och }} \quad \beta _{k}=m_{k+1}/\delta _{k}$ .

   Om antingen $\alpha _{k}$ eller $\beta _{k}$ är negativ, är indatapunkterna inte strikt monotona, och ${\ displaystyle (x_{k},\,y_{k})}$ är ett lokalt extremum. I sådana fall bitvisa monotona kurvor fortfarande genereras genom att välja $m_{k}=0\,$ om $\alpha _{k}<0$ eller $m_{k+1}=0\,$ om $\beta _{k}<0$ , även om strikt monotoni inte är möjlig globalt.
   
   
5. För att förhindra överskjutning och säkerställa monotonitet måste minst ett av följande tre villkor vara uppfyllt:

(a) funktionen

$\phi _{k}=\alpha _{k}- {\frac {(2\alpha _{k}+\beta _{k}-3)^{2}}{3(\alpha _{k}+\beta _{k}-2)}}>0 \,$ , eller

(b) ${\displaystyle \alpha _{k}+2\beta _{k}-3\leq 0\,} ,$ eller (

c) $2\alpha _{k}+\beta _{k}-3\leq 0\,$ .

Endast villkor (a) är tillräckligt för att säkerställa strikt monotoni: $\phi _{k}$ måste vara positivt.

Ett enkelt sätt att uppfylla denna begränsning är att begränsa vektorn $(\alpha _{k},\,\beta _{k})$ till en cirkel med radie 3. Det vill säga, om ${\displaystyle \alpha _{k}^{2}+\beta _{k}^{2}>9\,} ,$ ställ sedan in

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

och skala om tangenterna via

$m_{k}=\tau _{k}\,\alpha _{k}\ ,\delta _{k}\quad {\text{ och }}\quad m_{k+1}=\tau _{k}\,\beta _{k}\,\delta _{k}\,$ .

Alternativt är det tillräckligt att begränsa $\alpha _{k}\leq 3$ och $\beta _{k}\leq 3\,$ . För att åstadkomma detta om ${\displaystyle \alpha _{k}>3{\text{ eller }}\beta _{k}>3\,} ställ sedan$ in $m_{k}=3\,\delta _{k}\,$ .

Kubisk interpolation

Efter förbearbetningen ovan är utvärdering av den interpolerade spline ekvivalent med cubic Hermite spline , med hjälp av data $x_{k}$ , $y_{k}$ och $m_{ k}$ för $k=1,\,\dots \,n$ .

För att utvärdera vid $x$ , hitta indexet $k$ i sekvensen där ${\displaystyle x} ,$ ligger mellan $x_{k}$ och ${\ displaystyle x_{k+1}}$ , det vill säga: $x_{k}\leq x\leq x_{k+1}$ . Beräkna

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

då är det interpolerade värdet

${ \displaystyle f_{\text{interpolerad}}(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)}$

där $h_{ii}$ är basfunktionerna för den kubiska Hermite spline .

Exempel implementering

Följande JavaScript- implementering tar en datamängd och producerar en monoton kubisk spline-interpolantfunktion:

  
    
    

   0
     
        
    
    
           
       0      0  
        
        
        
           0
             
    
    
    
       
       0       
               
          
    
         
    
       0      
          
          
    

        
        
    
    
             
       0       
                      
          
          
          
    
    
       0
       0       
                
           0 
              0
          
                         
                    
        
    
      

        
        
      
        
    
    
          
       0       
           
           
           
                   
              
              
              
              
              
          
    
        
        

    
      
        
             
        
        
        
        
           0      
            
                
               
                      
                       
             
                
                
                    
                  
                
            
        
          0 

        
             
        
                        
        
                   
                              
        
        
            
    









   0    
   0    
   
   

    
    0       
        




        
   
   
    0       
       0  0
       
       0
       
                

 /*  * Monotone kubisk spline-interpolation  * Användningsexempel listat längst ner; detta är ett fullt fungerande paket.  Till   exempel * kan detta utföras antingen på webbplatser som  * https://www.programiz.com/javascript/online-compiler/  * eller med nodeJS.  */  function  DEBUG  (  s  )  {  /* Avkommentera följande för att möjliggöra utförlig utmatning av lösaren: */  //console.log(s);  }  var  outputCounter  =  ;  var  createInterpolant  =  funktion  (  xs  ,  ys  )  {  var  i  ,  length  =  xs  .  längd  ;  // Deal with length issues  if  (  length  !=  ys  .  length  )  {  throw  'Behöver ett lika stort antal xs och ys.'  ;  }  if  (  längd  ===  )  {  return  funktion  (  x  )  {  return  ;  };  }  if  (  längd  ===  1  )  {  // Impl: Förberäkning av resultatet förhindrar problem om ys muteras senare och tillåter sophämtning av ys //  Impl: Unary plus konverterar korrekt värden till tal  var  result  =  +  ys  [  ];  return  funktion  (  x  )  {  returnera  resultat  ;  };  }  // Ordna om xs och ys så att xs sorteras  var  indexes  =  [];  för  (  i  =  ;  i  <  längd  ;  i  ++  )  {  index  .  trycka  (  i  );  }  index  .  sort  (  funktion  (  a  ,  b  )  {  return  xs  [  a  ]  <  xs  [  b  ]  ?  -  1  :  1  ;  });  var  oldXs  =  xs  ,  oldYs  =  ys  ;  // Impl: Att skapa nya arrayer förhindrar också problem om inmatningsarrayerna muteras senare  xs  =  [];  ys  =  [];  // Impl: Unary plus konverterar korrekt värden till tal  för  (  i  =  ;  i  <  längd  ;  i  ++  )  {  xs  [  i  ]  =  +  oldXs  [  indexerar  [  i  ]];  ys  [  i  ]  =  +  oldYs  [  indexerar  [  i  ]];  }  DEBUG  (  "debug: xs = [ "  +  xs  +  " ]"  )  DEBUG  (  "debug: ys = [ "  +  ys  +  " ]"  )  // Få på varandra följande skillnader och lutningar  var  dys  =  [],  dxs  =  [] ,  ms  =  [];  för  (  i  =  ;  i  <  längd  -  1  ;  i  ++  )  {  var  dx  =  xs  [  i  +  1  ]  -  xs  [  i  ],  dy  =  ys  [  i  +  1  ]  -  ys  [  i  ];  dxs  [  i  ]  =  dx  ;  dys  [  i  ]  =  dy  ;  ms  [  i  ]  =  dy  /  dx  ;  }  // Få grad-1 koefficienter  var  c1s  =  [  ms  [  ]];  för  (  i  =  ;  i  <  dxs  .  längd  -  1  ;  i  ++  )  {  var  m  =  ms  [  i  ],  mNext  =  ms  [  i  +  1  ];  if  (  m  *  mNext  <=  )  {  c1s  [  i  ]  =  ;  }  else  {  var  dx_  =  dxs  [  i  ],  dxNext  =  dxs  [  i  +  1  ],  common  =  dx_  +  dxNext  ;  c1s  [  i  ]  =  3  *  common  /  ((  common  +  dxNext  )  /  m  +  (  common  +  dx_  )  /  mNext  );  }  }  c1s  .  push  (  ms  [  ms  .  längd  -  1  ]);  DEBUG  (  "debug: dxs = [ "  +  dxs  +  " ]"  )  DEBUG  (  "debug: ms = [ "  +  ms  +  " ]"  )  DEBUG  (  "debug: c1s.length = "  +  c1s  .  length  )  DEBUG  (  " debug: c1s = [ "  +  c1s  +  " ]"  )  // Få grad-2 och grad-3 koefficienter  var  c2s  =  [],  c3s  =  [];  för  (  i  =  ;  i  <  cls  .  längd  -  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  [  i  ]  =  (  m_  -  c1  -  common_  )  *  invDx  ;  c3s  [  i  ]  =  common_  *  invDx  *  invDx  ;  }  DEBUG  (  "debug: c2s = [ "  +  c2s  +  " ]"  )  DEBUG  (  "debug: c3s = [ "  +  c3s  +  " ]"  )  // Returnera interpolant funktion  returfunktion  (  x  )  {  //  Punkten längst till höger i datasetet ska ge ett exakt resultat  var  i  =  xs  .  längd  -  1  ;  // Avkommentera följande för att endast returnera det interpolerade värdet.  //if (x == xs[i]) { returnera ys[i]; }   // Sök efter intervallet x är i, returnera motsvarande y om x är ett av de ursprungliga xs  var  low  =  ,  mid  ,  high  =  c3s  .  längd  -  1  ;  while  (  låg  <=  hög  )  {  mid  =  Math  .  golv  (  0,5  *  (  låg  +  hög  ));  var  xHere  =  xs  [  miden  ];  if  (  xHere  <  x  )  {  low  =  mid  +  1  ;  }  else  if  (  xHere  >  x  )  {  high  =  mid  -  1  ;  }  else  {  // Avkommentera följande för att endast returnera det interpolerade värdet.  //retur ys[miden];  låg  =  c3s  .  längd  -  1  ;  hög  =  medel  ;  bryta  ;  }  }  i  =  Matte  .  max  (  ,  hög  );  // Interpolera  var  diff  =  x  -  xs  [  i  ];  outputCounter  ++  ;  var  interpolatedValue  =  ys  [  i  ]  +  diff  *  (  c1s  [  i  ]  +  diff  *  (  c2s  [  i  ]  +  diff  *  c3s  [  i  ]));  // Värdet på interpolatorns derivata vid denna punkt.  var  derivativeValue  =  c1s  [  i  ]  +  diff  *  (  2  *  c2s  [  i  ]  +  diff  *  3  *  c3s  [  i  ]);  DEBUG  (  "debug: #"  +  outputCounter  +  ". x = "  +  x  +  ". i = "  +  i  +  ", diff = "  +  diff  +  ", interpolatedValue = "  +  interpolatedValue  +  ", derivativeValue = "  +  derivativeValue  ) ;  // Avkommentera följande för att endast returnera det interpolerade värdet.  // returnera interpolatedValue;  return  [  interpoleratVärde  ,  derivatVärde  ];  };  };  /*  Användningsexempel nedan kommer att uppskatta x^2 för 0 <= x <= 4.  Kommandoradsanvändningsexempel (kräver installation av nodejs):  node monotone-cubic-spline.js  */  var  X  =  [  ,  1  ,  2  ,  3  4  ]  ;  var  F  =  [  ,  1  ,  4  ,  9  ,  16  ];  var  f  =  skapa Interpolant  (  X  ,  F  );  var  N  =  X  .  längd  ;  konsol  .  log (  "  # BLOCK 0 :: Data för monotone-cubic-spline.js" )  ;  konsol  .  log  (  "X"  +  "\t"  +  "F"  );  for  (  var  i  =  ;  i  <  N  ;  i  +=  1  )  {  konsol  .  log  (  F  [  i  ]  +  '\t'  +  X  [  i  ]);  }  konsolen  .  log  (  ""  );  konsol  .  log  (  ""  );  konsol  .  log  (  "# BLOCK 1 :: Interpolerad data för monotone-cubic-spline.js" )  ;  konsol  .  log  (  " x "  +  "\t\t"  +  " P(x) "  +  "\t\t"  +  "dP(x)/dx "  );  var  meddelande  =  ''  ;  var  M  =  25  ;  för  (  vari  =  ;  i  <  =  M  ;  i  +  =  1  )  {  varx  =  X  [  ]  +  (  X  [  N  -  1  ]  -X  [  ]  )  *  i  /  M  ;  var  rvals  =  f  (  x  );  var  P  =  rvals  [  ];  var  D  =  rvals  [  1  ];  meddelande  +=  x  .  toPrecision  (  15  )  +  '\t'  +  P  .  toPrecision  (  15  )  +  '\t'  +  D  .  toPrecision  (  15  )  +  '\n'  ;  }  konsolen  .  log  (  meddelande  ); 

• Fritsch, FN; Carlson, RE (1980). "Monotone bitvis kubisk interpolation". SIAM Journal on Numerical Analysis . SIAM. 17 (2): 238–246. doi : 10.1137/0717021 .
• Dougherty, RL; Edelman, A.; Hyman, JM (april 1989). "Positivitets-, monotoni- eller konvexitetsbevarande kubisk och kvintisk Hermite-interpolation" . Beräkningsmatematik . 52 (186): 471–494. doi : 10.2307/2008477 .

externa länkar

• GPLv 2 licensierad C++ implementering: MonotCubicInterpolator.cpp MonotCubicInterpolator.hpp