如何使用blitz++

Question

我是c++的初学者。我学习c++的重点是做科学的计算。我想使用blitz++库。我试图解决rk4方法，但我没有得到代码的内部工作原理(我知道rk4算法)。

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

这是我的问题？

1. 在rhs_eval x中，y只是值。但是dydx是指针。因此rhs_eval的输出值将被分配给y。不需要返回任何东西。我说得对吗？
2. int n = y.extent(0)是做什么的？在注释n中，它是耦合方程的数目。范围(0)的含义是什么？范围有多大？那是什么“0”？是第一个元素的大小吗？
3. 我如何打印'y‘的值？格式是什么？我想从rk4中获得y的值，方法是从main调用它。那就印出来。

我使用msvs2019编译了blitz++，用cmake使用这些指令- 说明

我从这里得到了密码- 只给出函数

Answer 1

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

另一件很棒的事情是，不用编译就可以使用blitz++。在Visual 2019中，展开{project name}，右键单击"references“和"Manage NuGet Packages”搜索blitz++。下载吧。没有额外的链接或其他必须做。

Answer 2

1. 是的，更改y也是通过引用传递的。指针使用*或指针模板，引用使用&。
2. 你的向量有一维或扩展。一般说来，Array<T,n>是n阶张量，n=2是一个矩阵。.extend(0)是一维的大小，具有基于零的索引.
3. 这是复杂的，没有很好的记录。我是说闪电战图书馆提供的设施。您可以手动打印组件。由于某些原因，如果第一个print命令被注释掉，我的版本会产生内存错误。

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