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