NEDISS: Network Diffusion and Synchronization Simulator
EulerSolver.h
1//
2// Created by m4zz31 on 2/11/21.
3//
4
5#ifndef CPPPROJCT_EULERSOLVER_H
6#define CPPPROJCT_EULERSOLVER_H
7
8#include "GeneralSolver.h"
9
10template <typename Equation>
11class EulerSolver {
12public:
13 EulerSolver(){};
14 void evolve(double t,
15 double h,
16 double a,
17 std::vector<double> &b,
18 std::vector<double> &c,
19 std::vector<double> &d,
20 Equation &E,
21 FlowSpecs &Specs,
22 std::vector<double> &RK1,
23 std::vector<double> &RK2,
24 std::vector<double> &RK3,
25 std::vector<double> &RK4,
26 double &answer,
27 double * P);
28 void Term1(double t,
29 double h,
30 double a,
31 std::vector<double> &b,
32 std::vector<double> &c,
33 std::vector<double> &d,
34 Equation &E,
35 FlowSpecs &Specs,
36 std::vector<double> &RK1,
37 double * P);
38 void Term2(double t,
39 double h,
40 double a,
41 std::vector<double> &b,
42 std::vector<double> &c,
43 std::vector<double> &d,
44 Equation &E,
45 FlowSpecs &Specs,
46 std::vector<double> &RK1,
47 std::vector<double> &RK2,
48 double * P);
49 void Term3(double t,
50 double h,
51 double a,
52 std::vector<double> &b,
53 std::vector<double> &c,
54 std::vector<double> &d,
55 Equation &E,
56 FlowSpecs &Specs,
57 std::vector<double> &RK1,
58 std::vector<double> &RK2,
59 std::vector<double> &RK3,
60 double * P);
61 void Term4(double t,
62 double h,
63 double a,
64 std::vector<double> &b,
65 std::vector<double> &c,
66 std::vector<double> &d,
67 Equation &E,
68 FlowSpecs &Specs,
69 std::vector<double> &RK1,
70 std::vector<double> &RK2,
71 std::vector<double> &RK3,
72 std::vector<double> &RK4,
73 double * P);
74};
75
76
77
78template <typename Equation>
79void EulerSolver<Equation>::evolve(double t,
80 double h,
81 double a,
82 std::vector<double> &b,
83 std::vector<double> &c,
84 std::vector<double> &d,
85 Equation &E,
86 FlowSpecs &Specs,
87 std::vector<double> &RK1,
88 std::vector<double> &RK2,
89 std::vector<double> &RK3,
90 std::vector<double> &RK4,
91 double &answer,
92 double * P){
93
94 // The Euler solver tries to use the deepest derivative as possible
95 if ((*(P+3)) != 0) {
96 answer = a + h * (RK1[0] +
97 h * (RK2[0] +
98 h * (RK3[0] +
99 h * RK4[0])));
100 } else if ((*(P+2)) != 0) {
101 answer = a + h * (RK1[0] +
102 h * (RK2[0] +
103 h * RK3[0]));
104 } else if ((*(P+1)) != 0) {
105 answer = a + h * (RK1[0] +
106 h * RK2[0]);
107 } else {
108 answer = a + h * RK1[0];
109 }
110}
111
112
113template <typename Equation>
114void EulerSolver<Equation>::Term1(double t,
115 double h,
116 double a,
117 std::vector<double> &b,
118 std::vector<double> &c,
119 std::vector<double> &d,
120 Equation &E,
121 FlowSpecs &Specs,
122 std::vector<double> &RK1,
123 double * P){
124 // Initialize auxiliary vectors to any value
125 std::vector<double> T1={0},T2={1},T3={0},T4={1};
126
127 // Values defined for Euler.
128 T1[0] = 0;
129 T2[0] = 1;
130 T3[0] = 0;
131 T4[0] = 1;
132 E.UpdateFlowSpecs(T1,T2,T3,T4,d.size());
133 E.Reset();
134 PRINTF_DBG("EulerSolver reports: before computing the field, Specs.result was %f\n",Specs.result);
135 E.Field(t,a,b,c,d);
136 PRINTF_DBG("EulerSolver reports: after computing the field, Specs.result was %f\n",Specs.result);
137 RK1[0] = Specs.result; // (*(P+0)) * h
138 PRINTF_DBG("EulerSolver reports: the written value to RK1[0] was %f\n",RK1[0]);
139
140
141}
142
143
144template <typename Equation>
145void EulerSolver<Equation>::Term2(double t,
146 double h,
147 double a,
148 std::vector<double> &b,
149 std::vector<double> &c,
150 std::vector<double> &d,
151 Equation &E,
152 FlowSpecs &Specs,
153 std::vector<double> &RK1,
154 std::vector<double> &RK2,
155 double * P){
156 // Initialize auxiliary vectors to any value
157 std::vector<double> T1={0},T2={1},T3={0},T4={1};
158
159 // Values defined for Euler.
160 // Here RK1 is the first derivative of the other variables :-)
161 T1[0] = 0;
162 T2[0] = 1;
163 T3[0] = 0;
164 T4[0] = 1;
165 E.UpdateFlowSpecs(T1,T2,T3,T4,d.size());
166 // Second Taylor Order
167 E.Reset();
168 E.d1Field(t,a,b,c,d);
169 RK2[0] = Specs.result; // (*(P+1)) * h * h / 2
170 PRINTF_DBG("a=%f,b=%f, b size is %d and the result was %f\n",a,b[0],b.size(), RK2[0]);std::cout<<std::flush;
171};
172
173template <typename Equation>
174void EulerSolver<Equation>::Term3(double t,
175 double h,
176 double a,
177 std::vector<double> &b,
178 std::vector<double> &c,
179 std::vector<double> &d,
180 Equation &E,
181 FlowSpecs &Specs,
182 std::vector<double> &RK1,
183 std::vector<double> &RK2,
184 std::vector<double> &RK3,
185 double * P){
186
187 // Initialize auxiliary vectors to any value
188 std::vector<double> T1={0},T2={1},T3={0},T4={1};
189
190 // Values defined for Euler.
191 // Here RK1 and RK2 are the first amd second derivatives of the other variables :-)
192 T1[0] = 0;
193 T2[0] = 1;
194 T3[0] = 0;
195 T4[0] = 1;
196 E.UpdateFlowSpecs(T1,T2,T3,T4,d.size());
197 // Third Taylor Order
198 E.Reset();
199 E.d2Field(t,a,b,c,d);
200 RK3[0] = Specs.result; // (*(P+2)) * h * h * h / 6
201};
202
203
204template <typename Equation>
205void EulerSolver<Equation>::Term4(double t,
206 double h,
207 double a,
208 std::vector<double> &b,
209 std::vector<double> &c,
210 std::vector<double> &d,
211 Equation &E,
212 FlowSpecs &Specs,
213 std::vector<double> &RK1,
214 std::vector<double> &RK2,
215 std::vector<double> &RK3,
216 std::vector<double> &RK4,
217 double * P){
218
219 // Initialize auxiliary vectors to any value
220 std::vector<double> T1={0},T2={1},T3={0},T4={1};
221
222 // Values defined for Euler.
223 // Here RK1, RK2 and RK3 are the first, second and third derivatives of the other variables :-)
224 T1[0] = 0;
225 T2[0] = 1;
226 T3[0] = 0;
227 T4[0] = 1;
228 E.UpdateFlowSpecs(T1,T2,T3,T4,d.size());
229 E.Reset();
230 E.d3Field(t,a,b,c,d);
231 RK4[0] = Specs.result; // (*(P+3)) * h * h * h * h / 24
232};
233
234#endif //CPPPROJCT_EULERSOLVER_H
EulerSolver
Definition: EulerSolver.h:11
FlowSpecs
Definition: GeneralDifferentialEquation.h:15
Generated by doxygen 1.9.2