C/C++ Trying to access constructor variables in class functions

[cppIStough]

EDIT: I don't know why the cpp code box below isn't appearing in proper format. I'll fix it if I can, but I don't see anyhting wrong with the syntax to display properly.

Okay, this was helpful! Here's what I have (I'll clean up the rounding errors, and to answer a few posts ago, I am familiar with convergence, stiff ODEs, and I actually have heard of JIT compilation, though only about a year ago; I'm learning C++ because I want to know it (already know Python fairly well, though I'm sure compared to many on here I'm a novice), which is why I am picking some easy math problems to get my hands dirty)

But here's what I have so far in 3 files though I'm getting errors which I'll list below the three files this program consists of:

rk4.h

#pragma once
#include <iostream>
#include <vector>

class RK4 {
private:
    double x0, y0, x_end, h;

public:
    std::vector<double> vec_x, vec_y;

    void solve();
};

rk4.cpp

#include "rk4.h"

RK4::RK4(double f_(double x, double y), double x0_, double y0_, double x_end_, double h_)
{
    f = f_;

    x0 = x0_;
    y0 = y0_;
    x_end = x_end_;
    h = h_;

    int n = (x_end - x0) / h;
    vec_x.resize(n);
    vec_y.resize(n);
    for (int i = 0; i < n; ++i)
    {
        vec_x[i ] = i * h;
        vec_y[i ] = y0;
    }
}

void RK4::solve()
{
    for (int i = 0; i < vec_x.size() - 1; ++i) {
        double k1 = f(vec_x[i ], vec_y[i ]);
        double k2 = f(vec_x[i ] + h / 2, vec_y[i ] + h * k1 / 2);
        double k3 = f(vec_x[i ] + h / 2, vec_y[i ] + h * k2 / 2);
        double k4 = f(vec_x[i ] + h, vec_y[i ] + h * k3);
        vec_y[i + 1] = vec_y[i ] + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
    }
}

ODE.cpp

#include "rk4.h"
#include <iostream>

double dydx_equals(double x, double y)
{
    return y;
}

int main() {

    double h = 0.01;
    double x0 = 0.;
    double y0 = 1.;
    double x_end = 1.;

    RK4 solver(dydx_equals, x0, y0, x_end, h);
    solver.solve();

    for (auto v : solver.vec_y)
    {
        std::cout << v << "\n";
    }
    return 0;
}

In rk4.cpp there is an error with the constructor stating "no instance of overloaded function RK4::RK4" which I'm confused on. Another error in this file is with the f = f_, which reads "identifier "f" is undefined, which is why I thought I may need to declare it, but am unsure how since the function is not a member of RK4. Any help would be awesome!

 

[pbuk]

EDIT: I don't know why the cpp code box below isn't appearing in proper format. I'll fix it if I can, but I don't see anyhting wrong with the syntax to display properly.

It's because of the : this is formatting code for italics here so it breaks the formatter. Your code has now been fixed up by a moderator but in future it is best to use j, i0, ii or whatever for index counts here.

Okay, this was helpful! Here's what I have (I'll clean up the rounding errors, and to answer a few posts ago, I am familiar with convergence, stiff ODEs, and I actually have heard of JIT compilation, though only about a year ago; I'm learning C++ because I want to know it (already know Python fairly well, though I'm sure compared to many on here I'm a novice), which is why I am picking some easy math problems to get my hands dirty)

That's all good. Another good way to get into C++ IMHO is by playing with a microcontroller kit - also adds experience in hardware and electronics to your portfolio (or gives you some fun playing with these things if they are already an interest).

In rk4.cpp there is an error with the constructor stating "no instance of overloaded function RK4::RK4" which I'm confused on.

That's because you haven't declared RK4:RK4 in your header file:
rk4.h

class RK4
{
public:
    std::vector<double> vec_x, vec_y;
    RK4 (double f_(double x, double y), double x0_, double y0_, double x_end_, double h_);
}

Another error in this file is with the f = f_, which reads "identifier "f" is undefined, which is why I thought I may need to declare it, but am unsure how since the function is not a member of RK4.

In your original code, f was a member of RK4, and the way you have written the constructor you are still trying to initialize it as a member so if you want to do it that way you should declare it in the header file.

But you don't have to do it that way, I was thinking of a pattern more like this:

// rk4.h
class RK4
{
private:
    // We don't need any of these private variables any more.
public:
    // We don't need to define the constructor any more because it doesn't need to do anything.
    std::vector<double> vec_x, vec_y;
    void solve (double f(double x, double y), double x0, double y0, double x_end, double h);
}

// rk4.cpp
    void solve(double f(double x, double y), double x0, double y0, double x_end, double h)
    {
        // Build the vectors.
        // Do the solution.
    }

 

[Mark44]

It's because of the : this is formatting code for italics here so it breaks the formatter. Your code has now been fixed up by a moderator but in future it is best to use j, i0, ii or whatever for index counts here.

It's been fixed now. An index variable i can be used, provided that you add a space before or after the letter i. That way the MathJax interpreter doesn't convert the 'i' in brackets to the BBCode delimiter for italics.

 

[cppIStough]

It's because of the : this is formatting code for italics here so it breaks the formatter. Your code has now been fixed up by a moderator but in future it is best to use j, i0, ii or whatever for index counts here.That's all good. Another good way to get into C++ IMHO is by playing with a microcontroller kit - also adds experience in hardware and electronics to your portfolio (or gives you some fun playing with these things if they are already an interest).That's because you haven't declared RK4:RK4 in your header file:
rk4.h

class RK4
{
public:
    std::vector<double> vec_x, vec_y;
    RK4 (double f_(double x, double y), double x0_, double y0_, double x_end_, double h_);
}

In your original code, f was a member of RK4, and the way you have written the constructor you are still trying to initialize it as a member so if you want to do it that way you should declare it in the header file.

But you don't have to do it that way, I was thinking of a pattern more like this:

// rk4.h
class RK4
{
private:
    // We don't need any of these private variables any more.
public:
    // We don't need to define the constructor any more because it doesn't need to do anything.
    std::vector<double> vec_x, vec_y;
    void solve (double f(double x, double y), double x0, double y0, double x_end, double h);
}

// rk4.cpp
    void solve(double f(double x, double y), double x0, double y0, double x_end, double h)
    {
        // Build the vectors.
        // Do the solution.
    }

Okay great, this makes a lot of sense. I'll polish up the roundoff errors but here's what we have (which works and makes a ton of sense: thanks!!!). Only difference, which I think you intended, was to have solve not output a void, but instead output the solution vector (or were you imagining this differently? Which is preferred to you?)

rk4.h

#pragma once
#include <iostream>
#include <vector>class RK4 {
public:
    std::vector<double> vec_x, vec_y;

    std::vector<double> solve(double f(double x, double y), double x0, double y0, double x_end, double h);
};

rk4.cpp

#include "rk4.h"std::vector<double> RK4::solve(double f(double x, double y), double x0, double y0, double x_end, double h)
{
    int n = (x_end - x0) / h;
    vec_x.resize(n);
    vec_y.resize(n);

    for (int ii = 0; ii < n; ++ii)
    {
        vec_x[ii] = ii * h;
        vec_y[ii] = y0;
    }

    for (int ii = 0; i < vec_x.size() - 1; ++ii) {
        double k1 = f(vec_x[ii], vec_y[ii]);
        double k2 = f(vec_x[ii] + h / 2, vec_y[ii] + h * k1 / 2);
        double k3 = f(vec_x[ii] + h / 2, vec_y[ii] + h * k2 / 2);
        double k4 = f(vec_x[ii] + h, vec_y[ii] + h * k3);
        vec_y[ii + 1] = vec_y[ii] + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
    }

    return vec_y;
}

ODE.cpp

#include "rk4.h"
#include <iostream>

double dydx_equals(double x, double y)
{
    return y;
}

int main() {

    double h = 0.01;
    double x0 = 0.;
    double y0 = 1.;
    double x_end = 1.;

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

    for (auto v : soln)
    {
        std::cout << v << "\n";
    }
    return 0;
}

So thinking about expanding on top of this, it seems it would make more sense to make the "grid" in a separate class from RK4 since it's independent of Runge Kutta. What do you think? I know this program is small and simple but I want to program it in an intelligent manner.

And to reiterate what you said above, passing a function name with no args into a constructor/function/etc is a pointer to the function?

 

[Mark44]

And to reiterate what you said above, passing a function name with no args into a constructor/function/etc is a pointer to the function?

Yes.

 

[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;
}