龙格库塔法圆周运动

Question

我被要求解决这个微分方程：

（x，y，vx，vy）

'=（vx，vy，vy，-vx）

它应该返回一个周期2*pi的圆周运动。我实现了该功能：

class FunzioneBase 
{
  public:
    virtual VettoreLineare Eval(double t, const VettoreLineare& v) const = 0; 
};
class Circonferenza: public FunzioneBase
{
  private:
    double _alpha;
  public:
    Circonferenza(double alpha) { _alpha = alpha; };
    void SetAlpha(double alpha) { _alpha = alpha; };
    virtual VettoreLineare Eval(double t, const VettoreLineare& v) const;
};
VettoreLineare Circonferenza::Eval(double t, const VettoreLineare& v) const
{
  VettoreLineare y(4);
  if (v.GetN() != 4) 
  {
    std::cout << "errore" << std::endl;
    return 0;
  };
  y.SetComponent(0, v.GetComponent(2));
  y.SetComponent(1, v.GetComponent(3));
  y.SetComponent(2, pow(pow(v.GetComponent(0), 2.) + pow(v.GetComponent(1), 2.), _alpha) * v.GetComponent(3));
  y.SetComponent(3, - pow(pow(v.GetComponent(0), 2.)  + pow(v.GetComponent(1), 2.), _alpha)) * v.GetComponent(2));
  return y;
};

其中_alpha等于0。现在，这适用于欧拉方法：如果我为2 * pi * 10积分这个常，给定初始条件(1, 0, 0, -1)，精度0.003，我得到的位置与1 ± 0.1范围内的(1, 0)相当，正如我们应该期望的那样。但是，如果我将相同的 ODE 与 Runge Kutta 方法（精度0.003，持续 2 * pi * 10 秒）集成，则实现如下：

class EqDifferenzialeBase 
{
  public:
    virtual VettoreLineare Passo (double t, VettoreLineare& x, double h, FunzioneBase* f) const = 0;
};
class Runge_Kutta: public EqDifferenzialeBase 
{
  public:
    virtual VettoreLineare Passo(double t, VettoreLineare& v, double h,  FunzioneBase* f) const;
};
VettoreLineare Runge_Kutta::Passo(double t, VettoreLineare& v, double h, FunzioneBase* _f) const
{ 
  VettoreLineare k1 = _f->Eval(t, v);
  VettoreLineare k2 = _f->Eval(t + h / 2., v + k1 *(h / 2.));
  VettoreLineare k3 = _f->Eval(t + h / 2., v + k2 * (h / 2.));
  VettoreLineare k4 = _f->Eval(t + h, v + k3 * h);
  VettoreLineare y = v + (k1 + k2 * 2. + k3 * 2. + k4) * (h / 6.);
  return y;
}

程序返回一个近似等于0.39的x位置，而对于 4 阶龙格库塔方法，理论上精度应该在 1E-6 左右。我检查了一下，发现Runge_Kutta的时期似乎几乎翻了两番（因为在2 * pi的时间间隔中，x从1到0.48），但我不明白为什么。这是我的主要内容：

VettoreLineare DatiIniziali (4);
Circonferenza* myCirc = new Circonferenza(0);
DatiIniziali.SetComponent(0, 1.);
DatiIniziali.SetComponent(1, 0.);
DatiIniziali.SetComponent(2, 0.);
DatiIniziali.SetComponent(3, -1.);
double passo = 0.003;
Runge_Kutta myKutta;
for(int i = 0; i <= 2. * M_PI / passo; i++)
{
  DatiIniziali = myKutta.Passo(0, DatiIniziali, passo, myCirc);
  cout << DatiIniziali.GetComponent(0) << endl;
};
cout << 1 - DatiIniziali.GetComponent(0) << endl;

提前谢谢你。

Answer 1

更新：发现一个错误：始终使用 -Wall 选项进行编译，以捕获编译器的所有警告和自动代码更正。然后你会发现

fff: In member function ‘virtual VettoreLineare Circonferenza::Eval(double, const VettoreLineare&) const’:
fff:xxx:114: error: invalid operands of types ‘void’ and ‘double’ to binary ‘operator*’
     y.SetComponent(3, - pow(pow(v.GetComponent(0), 2.)  + pow(v.GetComponent(1), 2.), _alpha)) * v.GetComponent(2));
                                                                                                                  ^

您在_alpha后过早关闭的地方，因此SetComponent void成为一个因素。

---

更新 II：发现第二个错误：在您的另一篇文章中给出了线性向量类的代码。在那里，与加法（operator+）相反，标量向量积（operator*(double)）正在修改调用实例。因此，在计算k2中，k1的分量乘以h/2。但是，这种修改后的k1（以及修改后的k2和k3）被用于组装结果y导致一些几乎完全无用的状态更新。

---

原始的快速原型：我可以告诉你，python中精简的裸体实现可以完美地工作

import numpy as np
def odeint(f,t0,y0,tf,h):
    '''Classical RK4 with fixed step size, modify h to fit
        the full interval'''
    N = np.ceil( (tf-t0)/h )
    h = (tf-t0)/N
    t = t0
    y = np.array(y0)
    for k in range(int(N)):
        k1 = h*np.array(f(t      ,y       ))
        k2 = h*np.array(f(t+0.5*h,y+0.5*k1))
        k3 = h*np.array(f(t+0.5*h,y+0.5*k2))
        k4 = h*np.array(f(t+    h,y+    k3))
        y = y + (k1+2*(k2+k3)+k4)/6
        t = t + h
    return t, y
def odefunc(t,y):
    x,y,vx,vy = y
    return [ vx,vx,vy,-vx ]
pi = 4*np.arctan(1);
print odeint(odefunc, 0, [1,0,0,-1], 2*pi, 0.003)

结尾为

(t, [ x,y,vx,vy]) = (6.2831853071794184,
    [  1.00000000e+00,  -6.76088739e-15,   4.23436503e-12,
    -1.00000000e+00])

不出所料。您将需要一个调试器或中间输出来查找计算出错的地方。