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