我如何使用 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
}

这是我的问题？

在 rhs_eval 中，x,y 只是值。但是 dydx 是指针。所以rhs_eval的输出值会赋值给y。不需要 return 任何东西。我说的对吗？
int n = y.extent(0)是做什么的？在评论中 n 表示它是耦合方程的数量。范围（0）的含义是什么。范围有什么作用？那个'0'是什么？是第一个元素的大小吗？
如何打印 'y' 的值？格式是什么？我想通过从 main 调用它来从 rk4 获取 y 的值。然后打印出来。

我使用 MSVS 2019 和 cmake 使用这些指令编译了 blitz++—— Instruction

我从这里得到代码- only the function is given

Answer 1

是的，也更改y以通过引用传递。指针带有 * 或指针模板，引用带有 &.
您的矢量有 1 个维度或扩展。一般来说，Array<T,n> 是 n 阶张量，n=2 是矩阵。 .extend(0)是第一个维度的大小，索引为zero-based。
这很复杂而且没有很好的记录。我指的是 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;
}

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

另一件好事是，无需任何编译即可使用 blitz++。在 Visual Studio 2019 年，展开 {project name} 然后右键单击“references”和“Manage NuGet Packages”搜索 blitz++。下载它。无需添加额外的链接或其他操作。

我如何使用 blitz++

How do I use blitz++

c++

numerical-methods

blitz++