c++numerical-methodsblitz++

How do I use blitz++


I am a beginner in c++. My focus of learning c++ is to do scientific computation. I want to use blitz++ library. I am trying to solve rk4 method but I am not getting the inner workings of the code(I know rk4 algorithm)

#include <blitz/array.h>
#include <iostream>
#include <stdlib.h>
#include <math.h>

using namespace blitz;
using namespace std;

# This will evaluate the slopes. say if dy/dx = y, rhs_eval will return y.  
void rhs_eval(double x, Array<double, 1> y, Array<double, 1>& dydx)
{
    dydx = y;
}

void rk4_fixed(double& x, Array<double, 1>& y, void (*rhs_eval)(double, Array<double, 1>, Array<double, 1>&), double h)
{
    // Array y assumed to be of extent n, where n is no. of coupled equations
    int n = y.extent(0);

    // Declare local arrays
    Array<double, 1> k1(n), k2(n), k3(n), k4(n), f(n), dydx(n);

    // Zeroth intermediate step
    (*rhs_eval) (x, y, dydx);
    for (int j = 0; j < n; j++)
    {
        k1(j) = h * dydx(j);
        f(j) = y(j) + k1(j) / 2.;
    }

    // First intermediate step
    (*rhs_eval) (x + h / 2., f, dydx);
    for (int j = 0; j < n; j++)
    {
        k2(j) = h * dydx(j);
        f(j) = y(j) + k2(j) / 2.;
    }

    // Second intermediate step
    (*rhs_eval) (x + h / 2., f, dydx);
    for (int j = 0; j < n; j++)
    {
        k3(j) = h * dydx(j);
        f(j) = y(j) + k3(j);
    }

    // Third intermediate step
    (*rhs_eval) (x + h, f, dydx);
    for (int j = 0; j < n; j++)
    {
        k4(j) = h * dydx(j);
    }

    // Actual step
    for (int j = 0; j < n; j++)
    {
        y(j) += k1(j) / 6. + k2(j) / 3. + k3(j) / 3. + k4(j) / 6.;
    }
    x += h;
    
    return; # goes back to function. evaluate y at x+h without returning anything
}

int main()
{
    cout << y <<endl; # this will not work. The scope of y is limited to rk4_fixed
}

Here are my questions?

  1. In rhs_eval x,y are just values. But dydx is pointer. So rhs_eval's output value will be assigned to y. No need to return anything. Am i correct?

  2. What does int n = y.extent(0) do? In comment n is saying it's the number of coupled equation. What is the meaning of extent(0). what does extent do? what is that '0'? Is it the size of first element?

  3. How do I print the value of 'y'? what is the format? I want to get the value of y from rk4 by calling it from main. then print it.

I compiled blitz++ using MSVS 2019 with cmake using these instruction-- Instruction

I got the code from here- only the function is given


Solution

  • #include <blitz/array.h>
    #include <iostream>
    #include <cstdlib>
    
    using namespace blitz;
    using namespace std;
    
    /* This will evaluate the slopes. say if dy/dx = y, rhs_eval will return y.  */
    const double sig = 10; const double rho = 28; const double bet = 8.0 / 3;
    
    void lorenz(double x, Array<double, 1> y, Array<double, 1> &dydx)
    {
        /* y vector = x,y,z in components */
        dydx(0) = sig * (y(1) - y(0));
        dydx(1) = rho * y(0) - y(1) - y(0) * y(2);
        dydx(2) = y(0) * y(1) - bet * y(2);
    }
    
    void rk4_fixed(double& x, Array<double, 1>& y, void (*rhs_eval)(double, Array<double, 1>, Array<double, 1> &), double h)
    {
        int n = y.extent(0);
    
        Array<double, 1> k1(n), k2(n), k3(n), k4(n), f(n), dydx(n);
    
        (*rhs_eval) (x, y, dydx);
        for (int j = 0; j < n; j++)
        {
            k1(j) = h * dydx(j);
            f(j) = y(j) + k1(j) / 2.0;
        }
    
        (*rhs_eval) (x + h / 2., f, dydx);
        for (int j = 0; j < n; j++)
        {
            k2(j) = h * dydx(j);
            f(j) = y(j) + k2(j) / 2.;
        }
    
        (*rhs_eval) (x + h / 2., f, dydx);
        for (int j = 0; j < n; j++)
        {
            k3(j) = h * dydx(j);
            f(j) = y(j) + k3(j);
        }
    
        (*rhs_eval) (x + h, f, dydx);
        for (int j = 0; j < n; j++)
        {
            k4(j) = h * dydx(j);
        }
    
        for (int j = 0; j < n; j++)
        {
            y(j) += k1(j) / 6. + k2(j) / 3. + k3(j) / 3. + k4(j) / 6.;
        }
        x += h;
    }
    
    int main()
    {
        Array<double, 1> y(3);
        y = 1, 1, 1;
    
        double x = 0, h = 0.05;
        Array<double, 1> dydx(3);
        dydx = 0, 0, 0;
    
        
        for (int i = 0; i < 10; i++)
        {
            rk4_fixed(x, y, &lorenz, h);
            cout << x << " ,";
            for (int j = 0; j < 3; j++)
            {
                cout << y(j)<<" ";
            }
            cout << endl;
        }
    
        return 0;
    }
    

    Another great thing is, it is possible to use blitz++ without any compilation. In Visual Studio 2019, expand {project name} than right click "references" and "Manage NuGet Packages" search for blitz++. download it. No added extra linking or others have to be done.