# Howto simpleMatrixLeastSquareFit

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

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

## 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;
}

// ************************************************************************* //
```

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

## 2 Description

I want to do a curve-fit of the function