VC++中的Runge Kutta方法

你的代码中有几个问题...我会指出我所看到的。

class Runge_Kutta
{
public:
    float x[100];
    float y[200];
    float h;         // step size
    float K_1, K_2, K_3, K_4;
    float compute();  // function to compute data point
    float data();   //function to generate data signal 
    Runge_Kutta();
    ~Runge_Kutaa();
};

到目前为止，好吧...

Runge_Kutta::Runge_Kutta()
{
    x[100] = 0;            // input signal

这是你的第一个问题。x数组是 100 个浮点数，您正在尝试设置 101 个浮点数。请记住，C++使用从零开始的索引，因此唯一有效的索引是 0 到 99（包括 0 到 99）。通过写入x[100]你在数组之外写入...可能进入y，但你当然不应该这样做。

更有可能的是，您希望对整个数组进行零初始化。有不同的方法可以创建/初始化一系列值（例如向量、动态分配），但正如你所做的那样，最简单的方法可能是：

    for (int i = 0; i < 100; ++i) { x[i] = 0.0f; }

您最好创建一个命名常量来替换"幻数"100。 继续使用您的代码：

    y[200] = 0;            // output accumulator

上述x类似的问题。同样，我怀疑您正在尝试将整个数组零初始化;替换为：

    for (int i = 0; i < 200; ++i) { y[i] = 0.0f; }

然后：

    h = 0.2;

这将起作用;您可能会收到有关双精度到浮点数转换的警告，因为 2.0 被认为是双精度数，但您正在存储浮点数。您可以通过追加f来强制浮点文本。

    h = 0.2f;

继续：

    K_1 = K_2 = K_3 = K_4 = 0;
}
~Runge_Kutta::Runga_Kutta(){}
float Runge_Kutta::data()
{
    int i = 0;
    for(i=0 ; i<100 ; i++)
    {
        x[i] = 5*sin((2*3.14*5*i)*0.01);  /*x[i] = A*sin(2*Pi*f*i*Ts)
    }                                       Ts-sampling period, A-amplitude,
    return x[i];                            f-frequency*/
}

也许您实际上并没有使用多行注释（例如/* 和 */），但您实际上在这里注释掉了一些代码......足以使其不应编译。如果要在代码右侧创建块注释，请确保使用每行注释（例如//）。

假设你真的做对了，你在这里有一个更大的问题。首先声明 i ，将其设置为零，然后将其用作循环索引。就其本身而言，这不是问题，尽管在for内循环

for (int i = 0; i < 100; ++i)

最大的问题是在循环后访问x[i]。 此时i为 100，和以前一样，x[i] 在 x 数组之外，该数组只能从 0 到 99 进行索引。你的意思是要x[0]返回吗？你想要哪个 x？返回该特定值。

float Runge_Kutta::compute()
{
    int j = 0;
    int i = 0;
    for(i = 0 ; i<100 ; i++)       // indexing through the input samples
    {
        for(j = 0 ; j<i+1 ; j++) // compute data points btw i and i+1
        {
            K_1 = h*x[i];
            K_2 = h*((x[i]+(h/2) + K_1/2);
            K_3 = h*((x[i]+(h/2) + K_2/2);
            K_4 = h* ((x[i] + h) + K_3);
        }
        y[i] = (1/6)*(K_1 + 2*K_2 + 2*K_3 + K_4);  //new data value
    }
    return y[i];
}

您可能在本节中混淆了循环指示。（但我对龙格-库塔不够熟悉，不知道。 i应该是x和y的索引吗？ j在哪里习惯，应该在哪里？

最后，返回y[i]似乎会遇到与上一个循环类似的问题......实际上，您只是将y从 0 初始化为 99，然后返回y[100]（因为该函数末尾的 i == 100）。

如果你对我提出的任何观点感到困惑，我会尽量澄清;只是问。

编辑：如前所述，（1/6）将计算（使用整数除法）为零。将一个或两个常量替换为浮点数以确保浮点除法...例如 1.0/6

因此，您希望在调用data()或compute()后访问所有数据...但是现在您返回的是单个float而不是整个数组。 以下是一些选项。

首先，您已x可公开访问。通常不推荐，但如果这是学习或实验，因为它是 public ，您可以直接访问该数据。

Runge_Kutta Interpolate;
Interpolate.data();
Interpolate.compute();
for (int i = 0; i < 100; ++i) {
    printf("Data point x[%d] = %fn", i, Interpolate.x[i]);
}

您可以更改data()和compute()以返回指向数组的指针，但通常也不建议这样做：

float* Runge_Kutta::data()
{
    // the bulk of the function...
    return &x[0];
}
// Use:
float* pData = Interpolate.data();
printf("First data signal: %fn", pData[0]);

另一种选择是提供索引数据访问器：

class Runge_Kutta {
public:
    // everything else you already have...
    float getDataSignal(size_t i) const {
        return x[i];
    }
    float getDataPoint(size_t i) const {
        return y[i];
    }
};

甚至另一种选择是使用C++标准库中的容器，例如 vector ，其中包含许多方便的操作作为库的一部分。

您不应通过引入其他约束（例如必须插值数据集以获取系统函数）来使情况复杂化。

RK::dydx(x,y) 
{
    return 5*sin(2*M_PI*5*x);
}

并使用 RK4 食谱的正确转录

RK::step(x,y,h) 
{
    double K1 = h*dydx(x,y);
    double K2 = h*dydx(x+h/2, y+K1/2);
    double K3 = h*dydx(x+h/2, y+K2/2);
    double K4 = h*dydx(x+h  , y+K3  );
    return y+(K1+2*(K2+K3)+K4)/6
}

然后你可以循环

for(i=0;i<100-1; i++) 
{
   y[i+1] = step(x,y[i],h);
   x +=h;
}