1

この微分方程式を解くように求められました。

(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匹敵することがわかります。しかし、同じ ODE を次のように実装されたルンゲ・クッタの方法 (精度で数秒) と統合すると、次のようになります。(1, 0)1 ± 0.10.0032 * 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;
}

プログラムは、4 次のルンゲ クッタ法で、理論的には精度が必要な場合に、およそ にx等しい位置を返します。確認したところ、Runge_Kutta の周期はほぼ 4 倍になっているように見えます (経過中に からに到達するため) が、その理由がわかりません。これは私のメインの内容です：0.391E-62 * pix10.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;

前もって感謝します。

4

1 に答える 1

0

更新: 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て、voidのSetComponentが要因になる場所。

更新 II: 2 番目のエラーが特定されました:あなたの別の投稿では、線形ベクトル クラスのコードが示されています。ここでは、加算( operator+) とは対照的に、スカラー ベクトル積( operator*(double)) が呼び出し元のインスタンスを変更しています。したがって、 getk2のコンポーネントを計算する際に、k1が乗算されh/2ます。しかし、この変更されたk1(および変更された and も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])

予想通り。計算が間違っている場所を見つけるには、デバッガーまたは中間出力が必要です。

于 2016-01-23T18:22:34.880 に答える