如何使用闪电战++

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”的值？格式是什么？我想通过从 main 调用 rk4 来获取 y 的值。然后打印它。

我使用 MSVS 2019 和 cmake 使用这些指令编译了 blitz++ - 说明

我从这里得到了代码 - 只给出了功能

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 Studio 2019 中，展开{项目名称}，然后右键单击“引用”并“管理 NuGet 包”搜索 blitz++。下载它。无需添加额外的链接或其他操作。

Answer 2

1. 是的，也将 y 更改为通过引用传递。指针使用*或指针模板，引用使用&。

2. 您的向量具有 1 维或扩展。一般来说，Array 是一个 n 阶张量，对于 n=2 来说是一个矩阵。 .extend(0) 是第一个维度的大小，索引从零开始。

3. 这很复杂并且没有很好的文档记录。我指的是 Blitz 图书馆提供的设施。您只需手动打印组件即可。由于某种原因，如果第一个打印命令被注释掉，我的版本会产生内存错误。

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