Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
fe_series.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 
18 #include <deal.II/base/config.h>
19 
20 #include <deal.II/fe/fe_series.h>
21 
22 #include <iostream>
23 
24 
25 DEAL_II_NAMESPACE_OPEN
26 
27 namespace FESeries
28 {
29  std::pair<double, double>
30  linear_regression(const std::vector<double> &x, const std::vector<double> &y)
31  {
32  FullMatrix<double> K(2, 2), invK(2, 2);
33  Vector<double> X(2), B(2);
34 
35  Assert(x.size() == y.size(),
36  ExcMessage("x and y are expected to have the same size"));
37 
38  Assert(x.size() >= 2,
40  "at least two points are required for linear regression fit"));
41 
42  double sum_1 = 0.0, sum_x = 0.0, sum_x2 = 0.0, sum_y = 0.0, sum_xy = 0.0;
43 
44  for (unsigned int i = 0; i < x.size(); i++)
45  {
46  sum_1 += 1.0;
47  sum_x += x[i];
48  sum_x2 += x[i] * x[i];
49  sum_y += y[i];
50  sum_xy += x[i] * y[i];
51  }
52 
53  K(0, 0) = sum_1;
54  K(0, 1) = sum_x;
55  K(1, 0) = sum_x;
56  K(1, 1) = sum_x2;
57 
58  B(0) = sum_y;
59  B(1) = sum_xy;
60 
61  invK.invert(K);
62  invK.vmult(X, B, false);
63 
64  return std::make_pair(X(1), X(0));
65  }
66 } // namespace FESeries
67 
68 DEAL_II_NAMESPACE_CLOSE
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
std::pair< double, double > linear_regression(const std::vector< double > &x, const std::vector< double > &y)
Definition: fe_series.cc:30
static ::ExceptionBase & ExcMessage(std::string arg1)
void invert(const FullMatrix< number2 > &M)
#define Assert(cond, exc)
Definition: exceptions.h:1407