我被要求解决这个微分方程:
(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
。现在,这适用于欧拉方法:如果我将此 ODE 积分为2 * pi * 10
,给定初始条件(1, 0, 0, -1)
,并具有一定的0.003
精度,我得到的位置与(1, 0)
的范围内相当1 ± 0.1
,正如我们所期望的那样。但是,如果我将相同的 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;
}
对于四阶龙格库塔方法,该程序返回一个近似x
等于的位置,当理论上精度应该是,大约。我检查并发现 Runge_Kutta 的时期似乎几乎是一式四份(因为在一段时间内,从到),但我不明白为什么。这是我的主要内容:0.39
1E-6
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;
先感谢您。