CODE — C++ ODE Solver

A C++ solver for ordinary differential equations (ODEs) using piecewise quadratic collocation whilst utilizing dual-number automated differentiation. The solver finds approximate solutions to first-order ODEs of the form: (Which is essentially all first order ODEs :])

df/dt = g(f(t), t)

It functions by breaking the domain into finite intervals and fitting a quadratic polynomial, utilizing Newton-Rhapson iteration (Jacobian via AD) to approximate the ODE at specific collocation points.

---

How It Works :)

Mathematics Formulation and Explanation (for those who want to learn)

1. Interval Decomposition

The solution domain [T0, Max] is split into steps of size H. On each arbitrary interval [x, x+H], the solution is approximated as a quadratic polynomial:

F(z) = az² + bz + c

where z ∈ [-1, 1] is a local coordinate obtained via the linear map:

z = 2(t - T0) / H - 1

The Chebyshev-style (foreshadowing ;) )rescaling to [-1, 1] improves numerical conditioning. Furthermore, I just might add Gegenbauer polynomials altogether!

2. Collocation Conditions

Differentiating the quadratic and enforcing the ODE at z = -1 and z = 1 gives two residual equations:

R1(a, b) = 2a + b  − (H/2) · g(a + b + c,  T0 + H) = 0
R2(a, b) = −2a + b − (H/2) · g(a − b + c,  T0)     = 0

The third degree of freedom is fixed by the initial condition F(-1) = F0, which gives c = F0 + b − a.

3. Newton-Rhapson iteration

The 2×2 Jacobian of (R1, R2) with respect to (a, b) is computed using dual numbers (Dual<float>), by seeding one variable at a time:

a.eps = 1.f; b.eps = 0.f;  // compute ∂R1/∂a, ∂R2/∂a
a.eps = 0.f; b.eps = 1.f;  // compute ∂R1/∂b, ∂R2/∂b

Cramer's rule then solves for the Newton update (Δa, Δb) in closed form. This iterates 10 times per interval (convergence check not yet implemented).

4. Curve Assembly

Each solved quadratic is stored in a SmoothCurveFromQuadratics object. Evaluation at any t selects the correct interval and applies the local quadratic, this class also transforms each point t into z space.

---

Key Components

Dual<T> — a dual number re + eps·ε where ε² = 0.

Quadratic<T> — stores coefficients (a, b, c) and evaluates az² + bz + c at a given z.

SmoothCurveFromQuadratics — a piecewise curve assembled from a sequence of quadratics, one per interval. Handles the t → z coordinate transform on lookup.

SolveODE(F0, T0, H, Max) — the current main solver. Returns a SmoothCurveFromQuadratics approximating the ODE solution over [T0, Max].

---

Demo

The included demo solves the logistic equation:

df/dt = f(1 − f),   f(0) = 0.5

whose exact solution is the sigmoid:

f(t) = 1 / (1 + e^(−t))

The solver runs with step size H = 0.2 over t ∈ [0, 5] and prints the approximate value, exact value, and relative percentage error at each point.

Sample output:

...
f(0.5)=0.622277, real value: 0.622459, error: 0.0292897%
f(0.6)=0.645448, real value: 0.645656, error: 0.0322998%
f(0.7)=0.667961, real value: 0.668188, error: 0.033927%
f(0.8)=0.689732, real value: 0.689974, error: 0.0350833%
...

---

Getting Started

Prerequisites

• C++17-compatible compiler (GCC, Clang, or MSVC)
• Visual Studio (optional, .sln included)

Build with Visual Studio

1. Clone the repository and open the .sln file.
2. Build and run.

Build with g++/clang++

git clone https://github.com/vernthedragon/CODE.git
cd CODE
g++ -std=c++17 -O2 -o ode_solver main.cpp
./ode_solver

---

Extending the Solver

The ODE to solve is defined as a single function (in the future I will make this an argument to pass):

DualDefault g(DualDefault ft, DualDefault t) {
    return ft * (DualDefault(1.f) - ft); // replace with your equation
}

Change this to any first-order ODE in the form f'(t) = g(f(t), t) and update the call to SolveODE with your initial condition, step size, and domain.

Planned improvements (noted in source):

• Convergence check instead of fixed 10-iteration Newton loop
• Use products instead of divisions in Jacobian solve for performance

---

Roadmap

• Arbitrary-order polynomial approximation — the solver currently fits a quadratic (degree 2) on each interval. This will be generalized to polynomials of any order.

• Chebyshev & Gegenbauer Polynomials — Support for Chebyshev and Gegenbauer polynomial bases. These offer superior approximation properties over uniform spacing and are the natural fit for the [-1, 1] interval rescaling already implemented.

• PDE support — the solver will eventually be extended to handle partial differential equations.

• ??? — Physics simulation support maybe...

---

License

GPL-3.0 — see LICENSE for details.

---

Author

vernthedragon — GitHub automated differentiation — GitHub mathematics — GitHub