Skip to content

Commit ff5735b

Browse files
WenyinWeiWenyinWei
committed
readjust the solvers hierarchy in include folder
1 parent 029faa6 commit ff5735b

38 files changed

+2005
-2369
lines changed

‎examples/advanced_integrators_usage.cpp‎

Lines changed: 41 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -4,14 +4,8 @@
44
#include <iomanip>
55
#include <chrono>
66

7-
// Include all integrators
8-
#include <solvers/rk4_solver.hpp>
9-
#include <solvers/rk23_solver.hpp>
10-
#include <solvers/rk45_solver.hpp>
11-
#include <solvers/dop853_solver.hpp>
12-
#include <solvers/radau_solver.hpp>
13-
#include <solvers/bdf_solver.hpp>
14-
#include <solvers/lsoda_solver.hpp>
7+
// Include modern diffeq library
8+
#include <diffeq.hpp>
159

1610
// Test systems
1711

@@ -88,37 +82,37 @@ void demonstrate_exponential_decay() {
8882
// Test all integrators
8983
{
9084
std::vector<double> y = {1.0};
91-
RK4Integrator<std::vector<double>> integrator(exponential_decay);
85+
diffeq::integrators::ode::RK4Integrator<std::vector<double>> integrator(exponential_decay);
9286
time_integrator(integrator, y, t_start, dt, t_end, "RK4");
9387
}
9488

9589
{
9690
std::vector<double> y = {1.0};
97-
RK23Integrator<std::vector<double>> integrator(exponential_decay, 1e-6, 1e-9);
91+
diffeq::integrators::ode::RK23Integrator<std::vector<double>> integrator(exponential_decay, 1e-6, 1e-9);
9892
time_integrator(integrator, y, t_start, dt, t_end, "RK23");
9993
}
10094

10195
{
10296
std::vector<double> y = {1.0};
103-
RK45Integrator<std::vector<double>> integrator(exponential_decay, 1e-8, 1e-12);
97+
diffeq::integrators::ode::RK45Integrator<std::vector<double>> integrator(exponential_decay, 1e-8, 1e-12);
10498
time_integrator(integrator, y, t_start, dt, t_end, "RK45");
10599
}
106100

107101
{
108102
std::vector<double> y = {1.0};
109-
// DOP853Integrator<std::vector<double>> integrator(exponential_decay, 1e-3, 1e-6);
110-
// time_integrator(integrator, y, t_start, dt, t_end, "DOP853");
103+
diffeq::integrators::ode::DOP853Integrator<std::vector<double>> integrator(exponential_decay, 1e-3, 1e-6);
104+
time_integrator(integrator, y, t_start, dt, t_end, "DOP853");
111105
}
112106

113107
{
114108
std::vector<double> y = {1.0};
115-
BDFIntegrator<std::vector<double>> integrator(exponential_decay, 1e-6, 1e-9);
109+
diffeq::integrators::ode::BDFIntegrator<std::vector<double>> integrator(exponential_decay, 1e-6, 1e-9);
116110
time_integrator(integrator, y, t_start, dt, t_end, "BDF");
117111
}
118112

119113
{
120114
std::vector<double> y = {1.0};
121-
LSODAIntegrator<std::vector<double>> integrator(exponential_decay, 1e-6, 1e-9);
115+
diffeq::integrators::ode::LSODA<std::vector<double>> integrator(exponential_decay, 1e-6, 1e-9);
122116
time_integrator(integrator, y, t_start, dt, t_end, "LSODA");
123117
}
124118
}
@@ -133,7 +127,7 @@ void demonstrate_van_der_pol() {
133127

134128
{
135129
std::vector<double> y = {1.0, 0.0};
136-
RK4Integrator<std::vector<double>> integrator(
130+
diffeq::integrators::ode::RK4Integrator<std::vector<double>> integrator(
137131
[&vdp](double t, const std::vector<double>& y, std::vector<double>& dydt) {
138132
vdp(t, y, dydt);
139133
});
@@ -142,25 +136,26 @@ void demonstrate_van_der_pol() {
142136

143137
{
144138
std::vector<double> y = {1.0, 0.0};
145-
RK45Integrator<std::vector<double>> integrator(
139+
diffeq::integrators::ode::RK45Integrator<std::vector<double>> integrator(
146140
[&vdp](double t, const std::vector<double>& y, std::vector<double>& dydt) {
147141
vdp(t, y, dydt);
148142
}, 1e-6, 1e-9);
149143
time_integrator(integrator, y, t_start, dt, t_end, "RK45");
150144
}
151145

152-
{
153-
std::vector<double> y = {1.0, 0.0};
154-
RadauIntegrator<std::vector<double>> integrator(
155-
[&vdp](double t, const std::vector<double>& y, std::vector<double>& dydt) {
156-
vdp(t, y, dydt);
157-
}, 1e-6, 1e-9);
158-
time_integrator(integrator, y, t_start, dt, t_end, "Radau");
159-
}
146+
// Note: RadauIntegrator not implemented in current hierarchy
147+
// {
148+
// std::vector<double> y = {1.0, 0.0};
149+
// RadauIntegrator<std::vector<double>> integrator(
150+
// [&vdp](double t, const std::vector<double>& y, std::vector<double>& dydt) {
151+
// vdp(t, y, dydt);
152+
// }, 1e-6, 1e-9);
153+
// time_integrator(integrator, y, t_start, dt, t_end, "Radau");
154+
// }
160155

161156
{
162157
std::vector<double> y = {1.0, 0.0};
163-
BDFIntegrator<std::vector<double>> integrator(
158+
diffeq::integrators::ode::BDFIntegrator<std::vector<double>> integrator(
164159
[&vdp](double t, const std::vector<double>& y, std::vector<double>& dydt) {
165160
vdp(t, y, dydt);
166161
}, 1e-6, 1e-9);
@@ -169,13 +164,13 @@ void demonstrate_van_der_pol() {
169164

170165
{
171166
std::vector<double> y = {1.0, 0.0};
172-
LSODAIntegrator<std::vector<double>> integrator(
167+
diffeq::integrators::ode::LSODA<std::vector<double>> integrator(
173168
[&vdp](double t, const std::vector<double>& y, std::vector<double>& dydt) {
174169
vdp(t, y, dydt);
175170
}, 1e-6, 1e-9);
176171
double timing = time_integrator(integrator, y, t_start, dt, t_end, "LSODA");
177172
std::cout << " Final method: " <<
178-
(integrator.get_current_method() == LSODAIntegrator<std::vector<double>>::MethodType::ADAMS ?
173+
(integrator.get_current_method() == diffeq::integrators::ode::LSODA<std::vector<double>>::MethodType::ADAMS ?
179174
"Adams (non-stiff)" : "BDF (stiff)") << std::endl;
180175
}
181176
}
@@ -190,25 +185,25 @@ void demonstrate_lorenz_system() {
190185

191186
{
192187
std::vector<double> y = {1.0, 1.0, 1.0};
193-
RK4Integrator<std::vector<double>> integrator(lorenz_system);
188+
diffeq::integrators::ode::RK4Integrator<std::vector<double>> integrator(lorenz_system);
194189
time_integrator(integrator, y, t_start, dt, t_end, "RK4");
195190
}
196191

197192
{
198193
std::vector<double> y = {1.0, 1.0, 1.0};
199-
RK45Integrator<std::vector<double>> integrator(lorenz_system, 1e-8, 1e-12);
194+
diffeq::integrators::ode::RK45Integrator<std::vector<double>> integrator(lorenz_system, 1e-8, 1e-12);
200195
time_integrator(integrator, y, t_start, dt, t_end, "RK45");
201196
}
202197

203198
{
204199
std::vector<double> y = {1.0, 1.0, 1.0};
205-
// DOP853Integrator<std::vector<double>> integrator(lorenz_system, 1e-3, 1e-6);
206-
// time_integrator(integrator, y, t_start, dt, t_end, "DOP853");
200+
diffeq::integrators::ode::DOP853Integrator<std::vector<double>> integrator(lorenz_system, 1e-3, 1e-6);
201+
time_integrator(integrator, y, t_start, dt, t_end, "DOP853");
207202
}
208203

209204
{
210205
std::vector<double> y = {1.0, 1.0, 1.0};
211-
LSODAIntegrator<std::vector<double>> integrator(lorenz_system, 1e-8, 1e-12);
206+
diffeq::integrators::ode::LSODA<std::vector<double>> integrator(lorenz_system, 1e-8, 1e-12);
212207
time_integrator(integrator, y, t_start, dt, t_end, "LSODA");
213208
}
214209
}
@@ -226,28 +221,29 @@ void demonstrate_stiff_robertson() {
226221
std::cout << "Using shortened time range (t=1.0) for demonstration purposes." << std::endl;
227222

228223
// Only test implicit methods for this stiff system
229-
try {
230-
std::vector<double> y = {1.0, 0.0, 0.0};
231-
RadauIntegrator<std::vector<double>> integrator(robertson_kinetics, 1e-6, 1e-9);
232-
time_integrator(integrator, y, t_start, dt, t_end, "Radau");
233-
} catch (const std::exception& e) {
234-
std::cout << std::setw(15) << "Radau" << ": Failed - " << e.what() << std::endl;
235-
}
224+
// Note: RadauIntegrator not implemented in current hierarchy
225+
// try {
226+
// std::vector<double> y = {1.0, 0.0, 0.0};
227+
// RadauIntegrator<std::vector<double>> integrator(robertson_kinetics, 1e-6, 1e-9);
228+
// time_integrator(integrator, y, t_start, dt, t_end, "Radau");
229+
// } catch (const std::exception& e) {
230+
// std::cout << std::setw(15) << "Radau" << ": Failed - " << e.what() << std::endl;
231+
// }
236232

237233
try {
238234
std::vector<double> y = {1.0, 0.0, 0.0};
239-
BDFIntegrator<std::vector<double>> integrator(robertson_kinetics, 1e-6, 1e-9);
235+
diffeq::integrators::ode::BDFIntegrator<std::vector<double>> integrator(robertson_kinetics, 1e-6, 1e-9);
240236
time_integrator(integrator, y, t_start, dt, t_end, "BDF");
241237
} catch (const std::exception& e) {
242238
std::cout << std::setw(15) << "BDF" << ": Failed - " << e.what() << std::endl;
243239
}
244240

245241
try {
246242
std::vector<double> y = {1.0, 0.0, 0.0};
247-
LSODAIntegrator<std::vector<double>> integrator(robertson_kinetics, 1e-6, 1e-9);
243+
diffeq::integrators::ode::LSODA<std::vector<double>> integrator(robertson_kinetics, 1e-6, 1e-9);
248244
double timing = time_integrator(integrator, y, t_start, dt, t_end, "LSODA");
249245
std::cout << " Final method: " <<
250-
(integrator.get_current_method() == LSODAIntegrator<std::vector<double>>::MethodType::ADAMS ?
246+
(integrator.get_current_method() == diffeq::integrators::ode::LSODA<std::vector<double>>::MethodType::ADAMS ?
251247
"Adams (non-stiff)" : "BDF (stiff)") << std::endl;
252248
} catch (const std::exception& e) {
253249
std::cout << std::setw(15) << "LSODA" << ": Failed - " << e.what() << std::endl;
@@ -269,7 +265,7 @@ void demonstrate_adaptive_features() {
269265

270266
for (auto [rtol, atol] : tolerances) {
271267
std::vector<double> y = {1.0};
272-
RK45Integrator<std::vector<double>> integrator(exponential_decay, rtol, atol);
268+
diffeq::integrators::ode::RK45Integrator<std::vector<double>> integrator(exponential_decay, rtol, atol);
273269

274270
std::cout << "Tolerances: rtol = " << rtol << ", atol = " << atol << std::endl;
275271
time_integrator(integrator, y, t_start, dt, t_end, "RK45");
@@ -296,7 +292,7 @@ int main() {
296292
std::cout << "✓ RK23: Adaptive 3rd order with embedded error estimation" << std::endl;
297293
std::cout << "✓ RK45: Adaptive 5th order Dormand-Prince (scipy default)" << std::endl;
298294
std::cout << "✓ DOP853: High-accuracy 8th order for demanding problems" << std::endl;
299-
std::cout << "✓ Radau: Implicit 5th order for stiff systems" << std::endl;
295+
std::cout << "• Radau: Implicit 5th order for stiff systems (not yet implemented)" << std::endl;
300296
std::cout << "✓ BDF: Variable-order implicit multistep for stiff systems" << std::endl;
301297
std::cout << "✓ LSODA: Automatic method switching (Adams ↔ BDF)" << std::endl;
302298

‎examples/quick_test.cpp‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@ int main() {
2121
try {
2222
std::cout << "RK4 (fixed step): ";
2323
auto y = y0;
24-
auto integrator = RK4Integrator<std::vector<double>>(exponential_decay);
24+
auto integrator = diffeq::integrators::ode::RK4Integrator<std::vector<double>>(exponential_decay);
2525
integrator.set_time(t_start);
2626
integrator.integrate(y, dt, t_end);
2727
std::cout << std::setprecision(6) << y[0] << std::endl;
@@ -32,7 +32,7 @@ int main() {
3232
try {
3333
std::cout << "RK23 (adaptive): ";
3434
auto y = y0;
35-
auto integrator = RK23Integrator<std::vector<double>>(exponential_decay, 1e-3, 1e-6);
35+
auto integrator = diffeq::integrators::ode::RK23Integrator<std::vector<double>>(exponential_decay, 1e-3, 1e-6);
3636
integrator.set_time(t_start);
3737
integrator.integrate(y, dt, t_end);
3838
std::cout << std::setprecision(6) << y[0] << std::endl;
@@ -43,7 +43,7 @@ int main() {
4343
try {
4444
std::cout << "RK45 (adaptive): ";
4545
auto y = y0;
46-
auto integrator = RK45Integrator<std::vector<double>>(exponential_decay, 1e-3, 1e-6);
46+
auto integrator = diffeq::integrators::ode::RK45Integrator<std::vector<double>>(exponential_decay, 1e-3, 1e-6);
4747
integrator.set_time(t_start);
4848
integrator.integrate(y, dt, t_end);
4949
std::cout << std::setprecision(6) << y[0] << std::endl;
@@ -54,7 +54,7 @@ int main() {
5454
try {
5555
std::cout << "DOP853 (high-acc): ";
5656
auto y = y0;
57-
auto integrator = DOP853Integrator<std::vector<double>>(exponential_decay, 1e-3, 1e-6);
57+
auto integrator = diffeq::integrators::ode::DOP853Integrator<std::vector<double>>(exponential_decay, 1e-3, 1e-6);
5858
integrator.set_time(t_start);
5959
integrator.integrate(y, dt, t_end);
6060
std::cout << std::setprecision(6) << y[0] << std::endl;
@@ -65,7 +65,7 @@ int main() {
6565
try {
6666
std::cout << "BDF (stiff): ";
6767
auto y = y0;
68-
auto integrator = BDFIntegrator<std::vector<double>>(exponential_decay, 1e-3, 1e-6);
68+
auto integrator = diffeq::integrators::ode::BDFIntegrator<std::vector<double>>(exponential_decay, 1e-3, 1e-6);
6969
integrator.set_time(t_start);
7070
integrator.integrate(y, dt, t_end);
7171
std::cout << std::setprecision(6) << y[0] << std::endl;
@@ -76,11 +76,11 @@ int main() {
7676
try {
7777
std::cout << "LSODA (automatic): ";
7878
auto y = y0;
79-
auto integrator = LSODAIntegrator<std::vector<double>>(exponential_decay, 1e-3, 1e-6);
79+
auto integrator = diffeq::integrators::ode::LSODA<std::vector<double>>(exponential_decay, 1e-3, 1e-6);
8080
integrator.set_time(t_start);
8181
integrator.integrate(y, dt, t_end);
8282
std::cout << std::setprecision(6) << y[0] << " (Method: " <<
83-
(integrator.get_current_method() == LSODAIntegrator<std::vector<double>>::MethodType::ADAMS ?
83+
(integrator.get_current_method() == diffeq::integrators::ode::LSODA<std::vector<double>>::MethodType::ADAMS ?
8484
"Adams)" : "BDF)") << std::endl;
8585
} catch (const std::exception& e) {
8686
std::cout << "Failed: " << e.what() << std::endl;

‎examples/rk4_integrator_usage.cpp‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@
33
#include <array>
44
#include <iomanip>
55
#include <memory>
6-
#include "solvers/rk4_solver.hpp"
6+
#include <diffeq.hpp>
77

88
// Example 1: Simple exponential decay
99
// dy/dt = -k*y, where k is the decay constant
@@ -50,7 +50,7 @@ int main() {
5050
std::cout << "\n1. Exponential Decay (dy/dt = -0.5*y)" << std::endl;
5151
std::cout << "-------------------------------------" << std::endl;
5252

53-
RK4Integrator<std::vector<double>, double> decay_integrator(exponential_decay);
53+
diffeq::integrators::ode::RK4Integrator<std::vector<double>, double> decay_integrator(exponential_decay);
5454
std::vector<double> decay_state = {2.0}; // y(0) = 2
5555

5656
std::cout << "Time\tValue" << std::endl;
@@ -65,7 +65,7 @@ int main() {
6565
std::cout << "\n2. Lorenz Attractor (first 20 steps)" << std::endl;
6666
std::cout << "------------------------------------" << std::endl;
6767

68-
RK4Integrator<std::vector<double>, double> lorenz_integrator(lorenz_system);
68+
diffeq::integrators::ode::RK4Integrator<std::vector<double>, double> lorenz_integrator(lorenz_system);
6969
std::vector<double> lorenz_state = {1.0, 1.0, 1.0}; // Initial conditions
7070

7171
std::cout << "Time\tX\t\tY\t\tZ" << std::endl;
@@ -85,7 +85,7 @@ int main() {
8585
std::cout << "\n3. Damped Harmonic Oscillator (float precision)" << std::endl;
8686
std::cout << "----------------------------------------------" << std::endl;
8787

88-
RK4Integrator<std::array<float, 2>, float> oscillator_integrator(damped_oscillator);
88+
diffeq::integrators::ode::RK4Integrator<std::array<float, 2>, float> oscillator_integrator(damped_oscillator);
8989
std::array<float, 2> oscillator_state = {1.0f, 0.0f}; // x(0) = 1, v(0) = 0
9090

9191
std::cout << "Time\tPosition\tVelocity" << std::endl;
@@ -104,7 +104,7 @@ int main() {
104104
std::cout << "\n4. Polymorphic Usage" << std::endl;
105105
std::cout << "-------------------" << std::endl;
106106

107-
auto integrator = std::make_unique<RK4Integrator<std::vector<double>, double>>(exponential_decay);
107+
auto integrator = std::make_unique<diffeq::integrators::ode::RK4Integrator<std::vector<double>, double>>(exponential_decay);
108108
AbstractIntegrator<std::vector<double>, double>* base_ptr = integrator.get();
109109

110110
std::vector<double> poly_state = {5.0};

‎examples/test_rk4_only.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@ int main() {
1111

1212
try {
1313
std::vector<double> y = {1.0};
14-
auto integrator = RK4Integrator<std::vector<double>>(exponential_decay);
14+
auto integrator = diffeq::integrators::ode::RK4Integrator<std::vector<double>>(exponential_decay);
1515
integrator.set_time(0.0);
1616
integrator.integrate(y, 0.1, 1.0);
1717
std::cout << "RK4 result: " << y[0] << " (expected: " << std::exp(-1.0) << ")" << std::endl;

‎examples/working_integrators_demo.cpp‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ int main() {
2020

2121
try {
2222
auto y = y0;
23-
RK4Integrator<std::vector<double>> integrator(exponential_decay);
23+
diffeq::integrators::ode::RK4Integrator<std::vector<double>> integrator(exponential_decay);
2424
integrator.set_time(t_start);
2525
integrator.integrate(y, dt, t_end);
2626
std::cout << "RK4 (4th order fixed): " << std::setprecision(6) << y[0] << std::endl;
@@ -30,7 +30,7 @@ int main() {
3030

3131
try {
3232
auto y = y0;
33-
RK23Integrator<std::vector<double>> integrator(exponential_decay, 1e-6, 1e-9);
33+
diffeq::integrators::ode::RK23Integrator<std::vector<double>> integrator(exponential_decay, 1e-6, 1e-9);
3434
integrator.set_time(t_start);
3535
integrator.integrate(y, dt, t_end);
3636
std::cout << "RK23 (3rd order adaptive): " << std::setprecision(6) << y[0] << std::endl;
@@ -40,7 +40,7 @@ int main() {
4040

4141
try {
4242
auto y = y0;
43-
RK45Integrator<std::vector<double>> integrator(exponential_decay, 1e-6, 1e-9);
43+
diffeq::integrators::ode::RK45Integrator<std::vector<double>> integrator(exponential_decay, 1e-6, 1e-9);
4444
integrator.set_time(t_start);
4545
integrator.integrate(y, dt, t_end);
4646
std::cout << "RK45 (5th order adaptive): " << std::setprecision(6) << y[0] << std::endl;
@@ -50,7 +50,7 @@ int main() {
5050

5151
try {
5252
auto y = y0;
53-
BDFIntegrator<std::vector<double>> integrator(exponential_decay, 1e-6, 1e-9);
53+
diffeq::integrators::ode::BDFIntegrator<std::vector<double>> integrator(exponential_decay, 1e-6, 1e-9);
5454
integrator.set_time(t_start);
5555
integrator.integrate(y, dt, t_end);
5656
std::cout << "BDF (stiff systems): " << std::setprecision(6) << y[0] << std::endl;
@@ -60,12 +60,12 @@ int main() {
6060

6161
try {
6262
auto y = y0;
63-
LSODAIntegrator<std::vector<double>> integrator(exponential_decay, 1e-6, 1e-9);
63+
diffeq::integrators::ode::LSODA<std::vector<double>> integrator(exponential_decay, 1e-6, 1e-9);
6464
integrator.set_time(t_start);
6565
integrator.integrate(y, dt, t_end);
6666
std::cout << "LSODA (automatic): " << std::setprecision(6) << y[0] <<
6767
" (" << (integrator.get_current_method() ==
68-
LSODAIntegrator<std::vector<double>>::MethodType::ADAMS ? "Adams" : "BDF") << ")" << std::endl;
68+
diffeq::integrators::ode::LSODA<std::vector<double>>::MethodType::ADAMS ? "Adams" : "BDF") << ")" << std::endl;
6969
} catch (const std::exception& e) {
7070
std::cout << "LSODA: Failed - " << e.what() << std::endl;
7171
}

0 commit comments

Comments
 (0)