Howto simpleMatrixLeastSquareFit

This is a very simple example of how to use the simpleMatrix class to solve a vector-matrix system using the LUsolve functionality

It does not require any dictionaries. I called it clduFoam and I execute it using clduFoam -e0 0 -e0 1.

The values used to fit the function are hardcoded because I wanted to keep it simple

1 Code

 
/*---------------------------------------------------------------------------*\
 Application
     clduFoam
 
 // base functions
 // 0 - 1
 // 1 - 1.0/T
 // 2 - log(T)
 // 3 - T^e
 
 Description
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "simpleMatrix.H"
 #include "OFstream.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Main program:
 
 int main(int argc, char *argv[])
 {
     argList::validOptions.insert("e0", "scalar");
     argList::validOptions.insert("e1", "scalar");
 
 #   include "setRootCase.H"
 
     // number of data-points
     label N = 4;
 
     // number of base functions
     label Nc = 4;
 
     // number of e-values between e0 and e1
     label Ne = 10;
 
     simpleMatrix<scalar> A(Nc);
     scalarField coeffs(Nc);
  
     scalar e0 = atof(args.options()["e0"].c_str());
     scalar e1 = atof(args.options()["e1"].c_str());
     scalar de = (e1-e0)/(Ne-1);
 
     scalarField T(N);
     scalarField pv(N);
     scalarField w(N);
 
     T[0] = 340.0; pv[0] = 0.133;   w[0] = 5.0;
     T[1] = 368.0; pv[1] = 1.33;    w[1] = 1.0;
     T[2] = 465.0; pv[2] = 1.0e+5;  w[2] = 100.0;
     T[3] = 700.0; pv[3] = 1.0e+8;  w[3] = 1.0e-0;
 
     scalarField logPv(log(pv));
     scalarField logT(log(T));
     scalarField Tinv(1.0/T);
 
     for(label it=0; it<Ne; it++)
     {
         scalar e = e0 + it*de;
         scalarField powT(pow(T,e));
    
         A.source()[0] = sum(logPv*w);
         A.source()[1] = sum(logPv*Tinv*w);
         A.source()[2] = sum(logPv*logT*w);
         A.source()[3] = sum(logPv*powT*w);
 
         A[0][0] = sum(w);
         A[0][1] = sum(Tinv*w);
         A[0][2] = sum(logT*w);
         A[0][3] = sum(powT*w);
    
         A[1][0] = sum(Tinv*w);
         A[1][1] = sum(Tinv*Tinv*w);
         A[1][2] = sum(logT*Tinv*w);
         A[1][3] = sum(powT*Tinv*w);
 
         A[2][0] = sum(logT*w);
         A[2][1] = sum(Tinv*logT*w);
         A[2][2] = sum(logT*logT*w);
         A[2][3] = sum(powT*logT*w);
 
         A[3][0] = sum(powT*w);
         A[3][1] = sum(Tinv*powT*w);
         A[3][2] = sum(logT*powT*w);
         A[3][3] = sum(powT*powT*w);
 
         coeffs = A.LUsolve();
         //Info << coeffs << endl;
         scalar err = 0.0;
         forAll(T,i)
         {
             scalar pc = coeffs[0] + coeffs[1]/T[i] + coeffs[2]*::log(T[i]) + coeffs[3]*::pow(T[i],e);
             err += ::pow(pc - logPv[i], 2);
             pc = ::exp(pc);
        
         }
         cout.precision(12);
         Info << "e = " << e << ", err = " << ::sqrt(err) << ",c = " << coeffs << endl;
     }
 
     forAll(T,i)
     {
         scalar pc = coeffs[0] + coeffs[1]/T[i] + coeffs[2]*::log(T[i]) + coeffs[3]*::pow(T[i],e1);
         pc = ::exp(pc);
         Info << "T = " << T[i] << ", pv = " << pv[i] << ", pvFit = " << pc << endl;
       
     }
 
     OFstream dataFile
     (
         "val.dat"
     );
  
     for(scalar Ti=300; Ti <= 710; Ti +=1)
     {
         scalar logpc = coeffs[0] + coeffs[1]/Ti + coeffs[2]*::log(Ti) + coeffs[3]*::pow(Ti,e1);
         scalar pc = ::exp(logpc);
         dataFile << Ti << " " << logpc << " " << pc << endl;
     }
 
     Info << "End\n" << endl;
     return 0;
 }
 
 
 // ************************************************************************* //
 

2 Description

I want to do a curve-fit of the function

$y(T) = \exp(c_0 + c_1/T + c_2 \log(T) + c_3 T^e)$

to a number of $P_v$ values.

I want to minimize the function

$S_e = \sum_i (y(T_i) - P_{v,i})^2$

with respect to the coefficients $c_i$.

But since the function is an exp, I will minimize the difference in log instead, i.e.

$S = \sum_{T_i} (\log(y(T_i)) - \log(P_{v,i}))^2$

Hence

$\frac{\partial S}{\partial c_i} = \frac{\partial}{\partial c_i}\sum_{T_i} [c_0 + c_1/T_i + c_2\log(T_i) + c_3T_i^e - log(P_{v,i})]^2 = 0$.

and we get

$\frac{\partial S}{\partial c_0} = 2\sum_{T_i} \{ [c_0 + c_1/T_i + c_2\log(T_i) + c_3T_i^e - log(P_{v,i})] \} = 0$

$\frac{\partial S}{\partial c_1} = 2\sum_{T_i} \{ [c_0 + c_1/T_i + c_2\log(T_i) + c_3T_i^e - log(P_{v,i})] (1/T_i)\} = 0$

$\frac{\partial S}{\partial c_2} = 2\sum_{T_i} \{ [c_0 + c_1/T_i + c_2\log(T_i) + c_3T_i^e - log(P_{v,i})] \log(T_i)\} = 0$

$\frac{\partial S}{\partial c_3} = 2\sum_{T_i} \{ [c_0 + c_1/T_i + c_2\log(T_i) + c_3T_i^e - log(P_{v,i})] T_i^e\} = 0$

which we can write as a vector matrix system $A c = b$

$(\sum_{T_i} 1)c_0 + (\sum_{T_i} 1/T_i)c_1 + (\sum_{T_i} \log(T_i))c_2 + (\sum_{T_i} T_i^e)c_3 = \sum_{T_i} log(P_{v,i})$

$(\sum_{T_i} 1/T_i)c_0 + (\sum_{T_i} (1/T_i)(1/T_i))c_1 + (\sum_{T_i} \log(T_i)/T_i)c_2 + (\sum_{T_i} T_i^e/T_i)c_3 = \sum_{T_i} log(P_{v,i})/T_i$

The matrix components should be clear from the code, where weights have also been added.

--Niklas 07:56, 21 January 2009 (CET)