nonlinear_least_squares.cpp

来自「算法的一些集合」· C++ 代码 · 共 33 行

#include "vs.h"
const double EPSILON = 1.e-12;
int main() {
   double weight[13] = {14.0/45.0, 64.0/45.0, 24.0/45.0, 64.0/45.0, 28.0/45.0,
                                   64.0/45.0, 24.0/45.0, 64.0/45.0, 28.0/45.0,
                                   64.0/45.0, 24.0/45.0, 64.0/45.0, 14.0/45.0};
   Quadrature qp(weight, 0.0, 1.0, 13);
   J d_l(1.0/12.0);      // per normalized length

   // B. Define Basis Functions
   H2 x((double*)0, qp),
   	phi = INTEGRABLE_VECTOR_OF_TANGENT_OF_TANGENT_BUNDLE("int, int, Quadrature",
            3/*vector size*/, 1/*spatial dim.*/, qp);
   phi[0] = 1.0-x.pow(2); phi[1] = 1.0-x.pow(4); phi[2] = 1.0-x.pow(6);
   H0 d2_phi = INTEGRABLE_VECTOR("int, Quadrature", 3, qp);
   for(int i = 0; i < 3; i++) d2_phi[i] = dd(phi)(i)[0][0]; // degenerated hessian matrix

   // C. Weighted-Residual Formulation 
   C0 c(3,(double*)0), delta_c(3,(double*)0);
   do { // Newton's Method
      H2 w = c * phi;
      H0	R = d(w)[0].pow(2)+(((H0)w)+sqrt(2.0))*dd(w)[0][0]-1.0,
         dR = 2.0*d(phi)(0)*d(w)[0]+((H0)phi)*dd(w)[0][0]+d2_phi*(((H0)w)+sqrt(2.0)),
         ddR = 2.0*(d(phi)(0)%d(phi)(0))+((H0)phi)%d2_phi+d2_phi%((H0)phi);
   	C0 df = ( dR*R ) | d_l,             // df / dc; gradient
   	   ddf = ( dR%dR+ddR*R ) | d_l;     // d^2f / dc^2; hessian
   	delta_c = -df /ddf; // Newton's formula
      c += delta_c;
      cout << c << ", " << norm(delta_c) << endl;                                // update increment
   } while((double)norm(delta_c ) > EPSILON); // check for convergence 
   cout << c << endl;
	return 0;
}