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
CodePudding user response:
Yes, change also
y
to be passed by reference. Pointer is with*
or a pointer template, reference is with&
.Your vector has 1 dimension or extend. In general
Array<T,n>
is a tensor of ordern
, forn=2
a matrix..extend(0)
is the size of the first dimension, with a zero-based index.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];
}
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;
}