我被要求解决这个微分方程:
((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
。现在,使用Euler的方法就可以很好地工作:如果我将2 * pi * 10
的ODE积分,给定初始条件(1, 0, 0, -1)
,且精度为0.003
,则可以得到在[...] (1, 0)
,正如我们所期望的。但是,如果我将同一ODE与Runge Kutta的方法集成在一起(精度为1 ± 0.1
,持续0.003
秒),则实现如下:
2 * pi * 10
如果理论上对于四阶Runge Kutta方法,精度应在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
位置,该位置近似等于0.39
。我检查发现,使用Runge_Kutta的周期似乎几乎翻了四倍(因为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;
选项进行编译,以捕获编译器的所有警告和自动代码更正。那么您会发现
-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
的void
成为一个因素。
Update II:识别出第二个错误:在SetComponent
中给出了线性向量类的代码。与addition(another post of yours)相反,scalar-vector product(operator+
)正在修改调用实例。因此,在计算operator*(double)
时,k2
的成分将与k1
相乘。但是,然后将此修改后的h/2
(以及修改后的k1
和k2
)用于组装结果k3
,从而导致几乎完全无用的状态更新。
原始快速原型:我可以告诉你,在python中精简的基本实现是完美的]]
y
以]结尾>
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)
如预期。您将需要调试器或中间输出来查找计算错误的地方。