How do I use blitz++

Question

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

Answer 1

1. Yes, change also y to be passed by reference. Pointer is with * or a pointer template, reference is with &.

2. Your vector has 1 dimension or extend. In general Array<T,n> is a tensor of order n, for n=2 a matrix. .extend(0) is the size of the first dimension, with a zero-based index.

3. This is complicated and not well documented. I mean the facilities provided by the Blitz library. You can just manually print the components. For some reason my version produces a memory error if the first print command is commented out.

#include <blitz/array.h>
#include <iostream>
#include <cstdlib>
//#include <cmath>

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];
*/
    /* use the comma operator */
    dydx = sig * (y[1] - y[0]), rho * y[0] - y[1] - y[0] * y[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)
{
    // 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);
    k1 = h * dydx; f=y+0.5*k1;
    

    // First intermediate step
    rhs_eval(x + 0.5*h, f, dydx);
    k2 = h * dydx; f =  y+0.5*k2;

    // Second intermediate step
    rhs_eval (x + 0.5*h, f, dydx);
    k3 = h * dydx; f=y+k3;
 
    // Third intermediate step
    rhs_eval (x + h, f, dydx);
    k4 = h * dydx;
 
    // Actual step
    y += k1 / 6. + k2 / 3. + k3 / 3. + k4 / 6.;
    x += h;
    
    return; //# goes back to function. evaluate y at x+h without returning anything
}

int main()
{
    Array<double, 1> y(3);
    y = 1,1,1;
    cout << y << endl;
    double x=0, h = 0.05;
    while(x<20) {
        rk4_fixed(x,y,lorenz,h);
        cout << x;
        for(int k =0; k<3; k++) {
            cout << ", "<< y(k);
        } 
        cout << endl;
    } 
    return 0;
}

Answer 2

#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.