Flow123d  release_2.2.0-914-gf1a3a4f
linear_ode_analytic.hh
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file linear_ode_analytic.hh
15  * @brief
16  */
17 
18 #ifndef LINEAR_ODE_ANALYTIC_H_
19 #define LINEAR_ODE_ANALYTIC_H_
20 
22 
23 /** @brief This class implements the analytic solution of a system of linear ODEs with constant matrix.
24  *
25  * The analytic solution can be obtained in the special case due to a physical nature of the problem.
26  * The problem is then solved only by a matrix multiplication.
27  *
28  * The assumption is made that the equations are independent. Each quantity is decreased (supposing
29  * negative diagonal) to \f$ e^{a_{ii} t} \f$. The decrement \f$ \left( 1-e^{a_{ii} t} \right) \f$
30  * is then distributed among other quantities according to the given fraction.
31  *
32  * In case of the decays and first order reactions the elements of the solution matrix are:
33  * \f{eqnarray*}{
34  * a_{ii} &=& e^{-\lambda_i t} \\
35  * a_{ji} &=& \left( 1-e^{-\lambda_i t} \right) b_{ji} \frac{M_j}{M_i}
36  * \f}
37  * where \f$ b_{ji} \f$ is the branching ratio of \f$ i \f$-th reactant and \f$ \frac{M_j}{M_i} \f$ is
38  * the fraction of molar masses.
39  *
40  * The fractions \f$ b_{ji} \frac{M_j}{M_i} \f$ are then obtained from the system matrix by dividing
41  * \f$ -\frac{a_{ji}}{a_{ii}} \f$.
42  *
43  * <B>Drawback:</B> These assumptions (equation independence) are adequate when very small time step is
44  * applied. This will lead to huge amount of evaluations of the exponential functions which can be expensive,
45  * so other numerical methods might be more appropriate.
46  * When the time step is large then the assumption is quite inadequate.
47  *
48  */
49 class LinearODEAnalytic : public LinearODESolver<LinearODEAnalytic>
50 {
51 public:
53 
54  /**
55  * Input record for class LinearODE_analytic.
56  */
57  static const Input::Type::Record & get_input_type();
58 
59  ///Default constructor is possible because the input record is not needed.
61 
62  /// Constructor from the input data.
64 
65  /// Destructor.
66  ~LinearODEAnalytic(void);
67 
68  void update_solution(arma::vec &init_vector, arma::vec &output_vec) override;
69 
70  bool evaluate_time_constraint(double &time_constraint) override;
71 
72 protected:
73  /**
74  * Computes the standard fundamental matrix.
75  */
76  void compute_matrix();
77 
78  /// The solution is computed only by a matrix multiplication (standard fundamental matrix).
79  arma::mat solution_matrix_;
80 
81 private:
82  /// Registrar of class to factory
83  static const int registrar;
84 
85 };
86 
87 #endif // LINEAR_ODE_ANALYTIC_H_
arma::mat solution_matrix_
The solution is computed only by a matrix multiplication (standard fundamental matrix).
This class implements the analytic solution of a system of linear ODEs with constant matrix...
LinearODESolverBase FactoryBaseType
Template class of the linear ODE solver.
static const Input::Type::Record & get_input_type()
bool evaluate_time_constraint(double &time_constraint) override
Estimate upper bound for time step. Return true if constraint was set.
Base class for linear ODE solver.
void update_solution(arma::vec &init_vector, arma::vec &output_vec) override
Updates solution of the ODEs system.
Accessor to the data with type Type::Record.
Definition: accessors.hh:292
static const int registrar
Registrar of class to factory.
~LinearODEAnalytic(void)
Destructor.
LinearODEAnalytic()
Default constructor is possible because the input record is not needed.
Record type proxy class.
Definition: type_record.hh:182