|
| 1 | + |
| 2 | +#include <fmu4cpp/fmu_base.hpp> |
| 3 | + |
| 4 | +#include <cmath> |
| 5 | + |
| 6 | +using namespace fmu4cpp; |
| 7 | + |
| 8 | +class VanDerPol : public fmu_base { |
| 9 | +public: |
| 10 | + FMU4CPP_CTOR(VanDerPol) { |
| 11 | + |
| 12 | + register_variable(real("x0", &x0_).setCausality(causality_t::OUTPUT).setInitial(initial_t::EXACT)); |
| 13 | + register_variable(real("x1", &x1_).setCausality(causality_t::OUTPUT).setInitial(initial_t::EXACT)); |
| 14 | + register_variable(real("mu", &mu_).setCausality(causality_t::PARAMETER).setVariability(variability_t::FIXED).setInitial(initial_t::EXACT)); |
| 15 | + |
| 16 | + register_variable(real("x0Der", &x0Der_).setCausality(causality_t::LOCAL).setInitial(initial_t::CALCULATED).setDerivative("x0", {"x0"})); |
| 17 | + register_variable(real("x1Der", &x1Der_).setCausality(causality_t::LOCAL).setInitial(initial_t::CALCULATED).setDerivative("x1", {"x0", "x1"})); |
| 18 | + |
| 19 | + VanDerPol::reset(); |
| 20 | + } |
| 21 | + |
| 22 | + bool do_step(double dt) override { |
| 23 | + if (dt <= 0.0) return false; |
| 24 | + |
| 25 | + // derivatives: |
| 26 | + // x' = y |
| 27 | + // y' = mu * (1 - x^2) * y - x |
| 28 | + auto dx = [&](double x, double y) { return y; }; |
| 29 | + auto dy = [&](double x, double y) { return mu_ * (1.0 - x * x) * y - x; }; |
| 30 | + |
| 31 | + // RK4 |
| 32 | + double k1x = dx(x0_, x1_); |
| 33 | + double k1y = dy(x0_, x1_); |
| 34 | + |
| 35 | + double x2 = x0_ + 0.5 * dt * k1x; |
| 36 | + double y2 = x1_ + 0.5 * dt * k1y; |
| 37 | + double k2x = dx(x2, y2); |
| 38 | + double k2y = dy(x2, y2); |
| 39 | + |
| 40 | + double x3 = x0_ + 0.5 * dt * k2x; |
| 41 | + double y3 = x1_ + 0.5 * dt * k2y; |
| 42 | + double k3x = dx(x3, y3); |
| 43 | + double k3y = dy(x3, y3); |
| 44 | + |
| 45 | + double x4 = x0_ + dt * k3x; |
| 46 | + double y4 = x1_ + dt * k3y; |
| 47 | + double k4x = dx(x4, y4); |
| 48 | + double k4y = dy(x4, y4); |
| 49 | + |
| 50 | + x0Der_ = k1x; |
| 51 | + x1Der_ = k1y; |
| 52 | + |
| 53 | + x0_ += (dt / 6.0) * (k1x + 2.0 * k2x + 2.0 * k3x + k4x); |
| 54 | + x1_ += (dt / 6.0) * (k1y + 2.0 * k2y + 2.0 * k3y + k4y); |
| 55 | + |
| 56 | + return true; |
| 57 | + } |
| 58 | + |
| 59 | + |
| 60 | + void reset() override { |
| 61 | + // default initial conditions |
| 62 | + x0_ = 2.0; |
| 63 | + x1_ = 0.0; |
| 64 | + mu_ = 1.0;// stiffness parameter; adjust if needed |
| 65 | + |
| 66 | + x0Der_ = 0.0; |
| 67 | + x1Der_ = 0.0; |
| 68 | + } |
| 69 | + |
| 70 | +private: |
| 71 | + double x0_{0.0}; |
| 72 | + double x1_{0.0}; |
| 73 | + double mu_{1.0}; |
| 74 | + |
| 75 | + double x0Der_{}; |
| 76 | + double x1Der_{}; |
| 77 | +}; |
| 78 | + |
| 79 | +model_info fmu4cpp::get_model_info() { |
| 80 | + model_info info; |
| 81 | + info.modelName = "VanDerPol"; |
| 82 | + info.description = "VanDerPol oscillator"; |
| 83 | + info.modelIdentifier = FMU4CPP_MODEL_IDENTIFIER; |
| 84 | + return info; |
| 85 | +} |
| 86 | + |
| 87 | +FMU4CPP_INSTANTIATE(VanDerPol); |
0 commit comments