CppNoddy  0.92
Loading...
Searching...
No Matches
Public Member Functions | List of all members
CppNoddy::ODE_IVP< _Type > Class Template Reference

A templated object for real/complex vector system of first-order ordinary differential equations. More...

#include <ODE_IVP.h>

Inheritance diagram for CppNoddy::ODE_IVP< _Type >:
CppNoddy::Uncopyable

Public Member Functions

 ODE_IVP (Equation< _Type > *equation_ptr, const double &x_init, const double &x_final, const std::size_t &num_of_steps)
 The class is defined by a vector function for the system. More...
 
 ~ODE_IVP ()
 
DenseVector< _Type > shoot4 (DenseVector< _Type > u)
 A fixed step 4th order Runge-Kutta method. More...
 
DenseVector< _Type > shoot45 (DenseVector< _Type > u, const double &tol, const double &h_init)
 A Runge-Kutta-Fehlberg integrator. More...
 
OneD_Node_Mesh< _Type > & get_mesh ()
 Return the history of the stepped solution. More...
 
unsigned & store_every ()
 Return a handle to the STORE_EVERY object. More...
 

Detailed Description

template<typename _Type>
class CppNoddy::ODE_IVP< _Type >

A templated object for real/complex vector system of first-order ordinary differential equations.

Definition at line 20 of file ODE_IVP.h.

Constructor & Destructor Documentation

◆ ODE_IVP()

template<typename _Type >
CppNoddy::ODE_IVP< _Type >::ODE_IVP ( Equation< _Type > *  equation_ptr,
const double &  x_init,
const double &  x_final,
const std::size_t &  num_of_steps 
)

The class is defined by a vector function for the system.

Parameters
equation_ptrA pointer to an inherited Equation object.
x_initThe starting point of the domain for the ODE.
x_finalThe end point of the domain for the ODE.
num_of_stepsA maximum/default number of steps to be taken.

Definition at line 19 of file ODE_IVP.cpp.

21 :
22 X_INIT(x1),
23 X_FINAL(x2),
24 H_INIT((x2 - x1) / num_of_points),
25 N(num_of_points),
26 p_EQUATION(ptr),
27 STORE_EVERY(1) {
28 p_EQUATION -> coord(0) = X_INIT;
29 }

◆ ~ODE_IVP()

template<typename _Type >
CppNoddy::ODE_IVP< _Type >::~ODE_IVP

Definition at line 32 of file ODE_IVP.cpp.

33 {}

Member Function Documentation

◆ get_mesh()

template<typename _Type >
OneD_Node_Mesh< _Type > & CppNoddy::ODE_IVP< _Type >::get_mesh

Return the history of the stepped solution.

Returns
A reference to the mesh solution

Definition at line 242 of file ODE_IVP.cpp.

242 {
243 return SOLN;
244 }

Referenced by main().

◆ shoot4()

template<typename _Type >
DenseVector< _Type > CppNoddy::ODE_IVP< _Type >::shoot4 ( DenseVector< _Type >  u)

A fixed step 4th order Runge-Kutta method.

Parameters
uAn NVector of initial values.
Returns
An NVector of the final values.

Definition at line 36 of file ODE_IVP.cpp.

36 {
37
38 double x = X_INIT;
39 const double h = H_INIT;
40 const double hby2 = h / 2.;
41 const double hby6 = h / 6.;
42 const double hby3 = h / 3.;
43 const int order = u.size();
44
45 DenseVector<_Type> z(order, 0.0), k1(order, 0.0), k2(order, 0.0), k3(order, 0.0), k4(order, 0.0);
46
47 DenseVector< double > coords;
48 std::vector< DenseVector<_Type> > values;
49
50 coords.push_back(x);
51 values.push_back(u);
52
53 for(unsigned i = 0; i < N; i++) {
54 // k1 = F(u,x)
55 p_EQUATION -> coord(0) = x;
56 p_EQUATION -> residual_fn(u, k1);
57 z = u + k1 * hby2;
58
59 x += hby2;
60
61 // k2 = F(z,xhh)
62 p_EQUATION -> coord(0) = x;
63 p_EQUATION -> residual_fn(z, k2);
64 z = u + k2 * hby2;
65
66 // k3 = F(z,xhh)
67 p_EQUATION -> coord(0) = x;
68 p_EQUATION -> residual_fn(z, k3);
69 z = u + k3 * h;
70
71 x += hby2;
72
73 // k4 = F(z,xh)
74 p_EQUATION -> coord(0) = x;
75 p_EQUATION -> residual_fn(z, k4);
76 u += k1 * hby6 + k2 * hby3 + k3 * hby3 + k4 * hby6;
77
78 if(i % STORE_EVERY == 0) {
79 coords.push_back(x);
80 values.push_back(u);
81 }
82
83 } //for loop stepping across domain
84
85 // construct the solution mesh stored in this object
86 SOLN = OneD_Node_Mesh<_Type>(coords, p_EQUATION -> get_order());
87 for(unsigned i = 0; i < coords.size(); ++i) {
88 // fill mesh
89 SOLN.set_nodes_vars(i, values[i]);
90 }
91 return u;
92 }
CppNoddy::DenseVector::push_back
void push_back(const _Type &fill)
A pass-thru definition of push_back.
Definition: DenseVector.h:310
CppNoddy::DenseVector::size
std::size_t size() const
A pass-thru definition to get the size of the vector.
Definition: DenseVector.h:330

◆ shoot45()

template<typename _Type >
DenseVector< _Type > CppNoddy::ODE_IVP< _Type >::shoot45 ( DenseVector< _Type >  u,
const double &  tol,
const double &  h_init 
)

A Runge-Kutta-Fehlberg integrator.

Parameters
uAn Nvector of initial values.
tolThe tolerance used in choosing the step length.
h_initThe initial step length.
Returns
An NVector of the final values.

Definition at line 95 of file ODE_IVP.cpp.

95 {
96 bool ok(false);
97 unsigned step = 0;
98 double x = X_INIT;
99 double h = h_init;
100 double c, diff;
101
102 static const double X2 = 1. / 4.;
103 static const double X3 = 3. / 8.;
104 static const double X4 = 12. / 13.;
105 static const double X5 = 1.;
106 static const double X6 = 1. / 2.;
107
108 static const double W21 = 1. / 4.;
109
110 static const double W31 = 3. / 32.;
111 static const double W32 = 9. / 32.;
112
113 static const double W41 = 1932. / 2197.;
114 static const double W42 = -7200. / 2197.;
115 static const double W43 = 7296. / 2197.;
116
117 static const double W51 = 439. / 216.;
118 static const double W52 = -8.;
119 static const double W53 = 3680. / 513.;
120 static const double W54 = -845. / 4104;
121
122 static const double W61 = -8. / 27.;
123 static const double W62 = 2.;
124 static const double W63 = -3544. / 2565.;
125 static const double W64 = 1859. / 4104.;
126 static const double W65 = -11. / 40.;
127
128 static const double U1 = 25. / 216.;
129 static const double U3 = 1408. / 2565.;
130 static const double U4 = 2197. / 4104.;
131 static const double U5 = -1. / 5.;
132
133 static const double Z1 = 16. / 135.;
134 static const double Z3 = 6656. / 12825.;
135 static const double Z4 = 28561. / 56430.;
136 static const double Z5 = -9. / 50.;
137 static const double Z6 = 2. / 55.;
138
139 const unsigned order = u.size();
140
141 DenseVector<_Type> z(order, 0.0), e(order, 0.0), k1(order, 0.0),
142 k2(order, 0.0), k3(order, 0.0), k4(order, 0.0), k5(order, 0.0), k6(order, 0.0);
143
144 DenseVector< double > coords;
145 std::vector< DenseVector<_Type> > values;
146
147 coords.push_back(x);
148 values.push_back(u);
149
150 do {
151 step += 1;
152 // k1 = F(u,x)
153 p_EQUATION -> coord(0) = x;
154 p_EQUATION -> residual_fn(u, k1);
155 k1 *= h;
156 z = u + k1 * W21;
157
158 // k2 = F(z,x+X2*h)
159 p_EQUATION -> coord(0) = x + X2 * h;
160 p_EQUATION -> residual_fn(z, k2);
161 k2 *= h;
162 z = u + k1 * W31 + k2 * W32;
163
164 // k3 = F(z,x+X3*h)
165 p_EQUATION -> coord(0) = x + X3 * h;
166 p_EQUATION -> residual_fn(z, k3);
167 k3 *= h;
168 z = u + k1 * W41 + k2 * W42 + k3 * W43;
169
170 // k4 = F(z,x+X4*h)
171 p_EQUATION -> coord(0) = x + X4 * h;
172 p_EQUATION -> residual_fn(z, k4);
173 k4 *= h;
174 z = u + k1 * W51 + k2 * W52 + k3 * W53 + k4 * W54;
175
176 // k5 = F(z,x+X5*h)
177 p_EQUATION -> coord(0) = x + X5 * h;
178 p_EQUATION -> residual_fn(z, k5);
179 k5 *= h;
180 z = u + k1 * W61 + k2 * W62 + k3 * W63 + k4 * W64 + k5 * W65;
181
182 // k6 = F(z,x+X6*h)
183 p_EQUATION -> coord(0) = x + X6 * h;
184 p_EQUATION -> residual_fn(z, k6);
185 k6 *= h;
186
187 e = k1 * U1 + k3 * U3 + k4 * U4 + k5 * U5;
188 z = k1 * Z1 + k3 * Z3 + k4 * Z4 + k5 * Z5 + k6 * Z6;
189
190 e -= z;
191
192 diff = e.inf_norm();
193
194 c = sqrt(sqrt(tol * h / (2 * diff)));
195 ok = true;
196
197 // is the first step ok? or does it need reducing?
198 if((step == 1) && (c < 1.0)) {
199 // step needs reducing so start from initial value again
200 ok = false;
201 step = 1;
202 }
203
204
205 if(ok) {
206 x += h;
207 u += z;
208
209 if(step % STORE_EVERY == 0) {
210 coords.push_back(x);
211 values.push_back(u);
212 }
213
214 }
215
216 h *= c;
217
218 if(x + h > X_FINAL) {
219 h = (X_FINAL - x);
220 }
221
222 if(step >= N) {
223 std::string problem;
224 problem = "The ODE.shoot45 method reached the maximum \n";
225 problem += "number of steps specified by the user. \n";
226 throw ExceptionRuntime(problem);
227 }
228
229 } while(std::abs(x - X_FINAL) > tol); // end loop stepping across domain
230
231 // construct the solution mesh stored in this object
232 SOLN = OneD_Node_Mesh<_Type>(coords, p_EQUATION -> get_order());
233 for(unsigned i = 0; i < coords.size(); ++i) {
234 // fill mesh
235 SOLN.set_nodes_vars(i, values[i]);
236 }
237
238 return u;
239 }

References h, CppNoddy::DenseVector< _Type >::push_back(), CppNoddy::DenseVector< _Type >::size(), and u.

Referenced by main().

◆ store_every()

template<typename _Type >
unsigned & CppNoddy::ODE_IVP< _Type >::store_every

Return a handle to the STORE_EVERY object.

Returns
A reference to the STORE_EVERY variable

Definition at line 247 of file ODE_IVP.cpp.

247 {
248 return STORE_EVERY;
249 }

Referenced by main().

The documentation for this class was generated from the following files:
Generated on Tue Jan 1 1980 00:00:00 for CppNoddy by   doxygen 1.9.5

© 2012

R.E. Hewitt