-
Notifications
You must be signed in to change notification settings - Fork 0
/
rhsOdeProblem.H
153 lines (113 loc) · 4.86 KB
/
rhsOdeProblem.H
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
# ifndef __RHS_ODE_PROBLEM_H__
# define __RHS_ODE_PROBLEM_H__
# include <string>
# include <functional>
# include <fstream>
namespace mg {
namespace numeric {
namespace ode {
/*-----------------------------------------------------------------------
* @brief Class RightHandSide - ODE problem,
*
* --> Type : single precision (float) , double precision (double)
*
* dy/dt = RHS
*
* @ Marco Ghiani Oct 2017 Glasgow UK
------------------------------------------------------------------------*/
template <typename Type = double>
class rhsOdeProblem {
//
//
//--
public:
rhsOdeProblem (const std::function<const Type(const Type,const Type)> numfun ,
const std::function<const Type(const Type,const Type)> exactfun ,
const Type,const Type,const Type,const Type,const std::string ) noexcept ;
rhsOdeProblem (const std::function<const Type(const Type,const Type)> numfun ,
const Type, const Type, const Type, const Type ) noexcept ;
virtual ~rhsOdeProblem() = default ;
rhsOdeProblem(const rhsOdeProblem &) = default ;
rhsOdeProblem(rhsOdeProblem&& ) = default ;
rhsOdeProblem& operator=(const rhsOdeProblem&) = default ;
rhsOdeProblem& operator=(rhsOdeProblem&& ) = default ;
std::function<const Type(const Type,const Type)> numericalFunction ;
std::function<const Type(const Type,const Type)> analiticalFunction ;
Type f(Type t, Type u) const noexcept { return numericalFunction(t,u); }
Type dfdt(Type t, Type u) const noexcept { return (f(t,u+eps)-f(t,u))/eps ; }
auto setRhs (std::function<Type(Type,Type)> numfun) noexcept { numericalFunction = numfun; }
auto setExact(std::function<Type(Type,Type)> exactfun) noexcept {analiticalFunction = exactfun;}
auto solveExact() noexcept ;
const Type t0() const noexcept { return _t0 ;}
const Type tf() const noexcept { return _tf ;}
const Type dt() const noexcept { return _dt ;}
const Type u0() const noexcept { return _u0 ;}
const std::string fname () const noexcept { return filename ;}
//---
private:
Type _t0 ;
Type _tf ;
Type _dt ;
Type _u0 ;
std::string filename ;
constexpr static Type eps = 1e-12 ;
};
/*
* Implementation
*/
template <typename Type>
rhsOdeProblem<Type>::rhsOdeProblem ( const std::function<const Type(const Type,const Type)> numfun ,
const std::function<const Type(const Type,const Type)> exactfun ,
const Type Ti,const Type Tf,const Type Dt,const Type U0,
const std::string fname
)
noexcept : numericalFunction{numfun} ,
analiticalFunction{exactfun} ,
_t0{Ti} ,
_tf{Tf} ,
_dt{Dt} ,
_u0{U0} ,
filename{fname}
{
solveExact();
}
template<typename Type>
rhsOdeProblem<Type>::rhsOdeProblem ( const std::function<const Type(const Type,const Type)> numfun ,
const Type Ti,const Type Tf,const Type Dt,const Type U0
)
noexcept : numericalFunction{numfun} ,
_t0{Ti} ,
_tf{Tf} ,
_dt{Dt} ,
_u0{U0}
{}
//- if exist ( and gives ) compute the
// numerical-exact solution
//
template<typename Type>
auto rhsOdeProblem<Type>::solveExact() noexcept {
const Type Ns = ( _tf -_t0 )/ _dt ;
std::ofstream fn( filename , std::ios::out);
if(!fn)
{
std::cerr << "Error opening file " << filename << " in rhsOdeProblem::solveExact()" << std::endl;
exit(-1);
}
else
{
//std::cout << Ns << std::endl ;
Type time = _t0 , yt = _u0 ;
fn << time << ' ' << yt << std::endl ;
for(std::size_t i=0 ; i < Ns ; i++)
{
time += _dt ;
yt = analiticalFunction(time,yt);
fn << time << ' ' << yt << std::endl ;
}
fn.close();
}
}
}//ode
}//numeric
}//mg
# endif