6 (応用例)観測点の多項式近似

6.1 概要

ここからはCPPLapackの機能の応用例を示していきます．最初の例は観測点の多項式近似です．一般に最小二乗法の方が通りがよいかもしれません．

6.2 アルゴリズムの解説

m個の観測点$(x_1,y_1)$〜$(x_n,y_n)$をm次の多項式近似する場合，

$$\begin{bmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^n\ 1 & x_2 & x... ... \left( \begin{array}{c} y_1\ y_2\ \vdots\ y_m \end{array} \right)$$

というAx=y型の方程式を解けば，多項式の係数$k_0〜k_n$を求められます． より詳しい解説は他を御覧下さい．

6.3 例題

4つの観測点 (0,-4.04), (1,-1.98), (2,2.02), (3,13.86) を3次の多項式で近似する．


コードを以下に示します．

//====================================================================[include]
#include ``cpplapack.h''

//======================================================================[solve]
void solve(int N,
           CPPL::dcovector& x, CPPL::dcovector& y)
{
  //// make dgematrix A ////
  CPPL::dgematrix A(x.l,N+1);
  for(int i=0; i<A.m; i++){for(int j=0; j<A.n; j++){
      A(i,j) = std::pow(x(i), (double)j);
  }}

  //// solve Ak=y ////
  A.dgels(y);
}

//=======================================================================[main]
/*! main */
int main(int argc, char** argv)
{
  //// make dcovector x,y ////
  CPPL::dcovector x(4), y(4);
  x(0)=0; y(0)=-4.04;
  x(1)=1; y(1)=-1.98;
  x(2)=2; y(2)=2.02;
  x(3)=3; y(3)=13.86;

  //// solve ////
  solve(3,x,y);

  //// print ////
  for(int i=0; i<3; i++){
    std::cout << ``k'' << i << ``='' << y(i) << std::endl;
  }

  return 0;
}

実行結果は以下の様になります．

k0=-4.04
k1=3.05667
k2=-1.98