Trying to access constructor variables in class functions

[cppIStough]

I can open a new thread if needed, but I thought since this question is right in line with what we spoke about above, perhaps it belongs here. I have the files

rk4.h

#pragma once
#include <vector>class ODESolverMethod {
public:
    double x0, y0, x_end;
    int n;
    explicit ODESolverMethod(double x0_, double y0_, double x_end_, int n_) {
        double x0 = x0_; double y0 = y0_; double x_end = x_end_; int n = n_;
    }

    virtual std::vector<std::pair<double, double>> 
        solve(double f(double x, double y)) = 0;

    std::vector<std::pair<double, double>> soln;
};

class RK4 : public ODESolverMethod {
public:
    explicit RK4(double x0_, double y0_, double x_end_, int n_) : ODESolverMethod(x0_, y0_, x_end_, n_) {
        double x0 = x0_; double y0 = y0_; double x_end = x_end_; int n = n_;
    }

    std::vector<std::pair<double, double>> solve(double f(double x, double y)) {
        double h = (x_end - x0) / n;

        soln.push_back(std::make_pair(x0, y0));

        auto y_i = y0;
        for (int i = 1; i < n + 1; ++i) {
            auto x_i = i * h;
            auto k1 = f(x_i, y_i);
            auto k2 = f(x_i + h / 2, y_i + h * k1 / 2);
            auto k3 = f(x_i + h / 2, y_i + h * k2 / 2);
            auto k4 = f(x_i + h, y_i + h * k3);
            y_i = y_i + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);

            soln.push_back(std::make_pair(x_i, y_i));
        }

        return soln;
    };

    std::vector<std::pair<double, double>> soln;
};

and ODE.cpp

#include "rk4.h"
#include <iostream>double dydx_equals(double x, double y)
{
    return x*x*x*exp(-2*x) - 2*y;
}

int main() {

    int n = 10;
    double x0 = 0.;
    double y0 = 1.;
    double x_end = 1.;

    RK4 solver(x0, y0, x_end, n);
    auto soln = solver.solve(dydx_equals);

    for (auto iter = soln.begin(); iter != soln.end(); ++iter)
    {
        std::cout << iter->first << ", " << iter->second << "\n";
    }
    return 0;
}

I'm trying to make a more general ODE library that can solve ODEs differently, say rk4 vs rk2 for starters. I think the best idea is to have a general class for the ODE method and the inherit from this. However, inheriting the constructor is a challenge. I've googled a ton and no matter what I do I can't get an output now that is not (0,0). If you see what I'm doing please help. I've stepped through code but can't see the issue.

 

[Filip Larsen]

You never set any of your class member because your constructor is only setting local stack variables you have named identical to your member variables (in line 10 in your header), so on a good day (e.g. during debug) n will have value zero when you solve. Similar, your line 22 is effectively doing nothing and can be removed.

Later: If I clean up your version I get something like

#include <cmath>
#include <iostream>
#include <vector>

class ODESolverMethod {
public:
    double x0, y0, x_end;
    int n;
    
    ODESolverMethod(double x0, double y0, double x_end, int n) 
        : x0(x0), y0(y0), x_end(x_end), n(n)
    {}

    virtual std::vector<std::pair<double, double>> 
        solve(double f(double x, double y)) = 0;
};

class RK4 : public ODESolverMethod {
public:
    RK4(double x0, double y0, double x_end, int n) 
        : ODESolverMethod(x0, y0, x_end, n) 
    {}

    std::vector<std::pair<double, double>> solve(double f(double x, double y)) {
        double h = (x_end - x0) / n;

        std::vector<std::pair<double, double>> soln;
        soln.push_back(std::make_pair(x0, y0));

        auto y_i = y0;
        std::cout << "n = " << n << std::endl;
        for (int i = 1; i < n + 1; ++i) {
            auto x_i = i * h;
            auto k1 = f(x_i, y_i);
            auto k2 = f(x_i + h / 2, y_i + h * k1 / 2);
            auto k3 = f(x_i + h / 2, y_i + h * k2 / 2);
            auto k4 = f(x_i + h, y_i + h * k3);
            y_i = y_i + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
            std::cout << x_i << ", " << y_i << std::endl;
            soln.push_back(std::make_pair(x_i, y_i));
        }

        return soln;
    }
};double dydx_equals(double x, double y)
{
    return x*x*x*std::exp(-2*x) - 2*y;
}

int main() {

    int n = 10;
    double x0 = 0.;
    double y0 = 1.;
    double x_end = 1.;

    RK4 solver(x0, y0, x_end, n);
    auto soln = solver.solve(dydx_equals);

    for (auto iter = soln.begin(); iter != soln.end(); ++iter)
    {
        std::cout << iter->first << ", " << iter->second << "\n";
    }
    return 0;
}

I have made a quick single-file version in a online compiler using C++20 for you to compare with, retaining a hierarchy (which is not really useful until yet another solver gets implemented) but here using lambda's (functional) instead of a function pointer and moving the initial state parameters to the solve call itself:

#include <cassert>
#include <cmath>
#include <iostream>
#include <functional>
#include <vector>
#include <tuple>

class Solver
{
    public:
        using ScalarSolution = std::vector<std::pair<double, double>>;
        using ScalarField = std::function<double(double, double)>;
    
        virtual ~Solver() = default;
    
        virtual ScalarSolution solve(double x0, double x1, double y0, double h, const ScalarField& f) const = 0;
};

class RK4Solver : public Solver
{
    public:
        ScalarSolution solve(double x0, double x1, double y0, double h, const ScalarField& f) const override
        {
            assert(h > 0);
            const auto n = static_cast<size_t>(std::ceil(std::abs(x1 - x0) / h));
            h = (x1 - x0) / n;
            auto y_i = y0;
        
            ScalarSolution sol;
            sol.reserve(n);
            sol.emplace_back(x0, y0);
            for (size_t i = 1u; i <= n; ++i)
            {
                const auto x_i = x0 + i * h;
                auto k1 = f(x_i, y_i);
                auto k2 = f(x_i + h / 2, y_i + h * k1 / 2);
                auto k3 = f(x_i + h / 2, y_i + h * k2 / 2);
                auto k4 = f(x_i + h, y_i + h * k3);
                y_i = y_i + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
                sol.emplace_back(x_i, y_i);
            }
            return sol;
        }
};int main()
{
    RK4Solver s;
    const auto sol = s.solve(0, 1, 1, 0.1, [](double x, double y) { return x*x*x*std::exp(-2*x) - 2*y; });
    for (const auto [x, y] : sol)
    {
        std::cout << x << ", " << y << "\n";
    }
    return 0;
}