Метод Рунге Кутты круговое движение

Question

Метод Рунге Кутты круговое движение

Меня попросили решить это дифференциальное уравнение:

(Х, у, 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 с методом Рунге Кутты (с 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 позиция, которая равна 0.39 Приблизительно, когда теоретически точность должна быть, для метода Рунге Кутты 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;

Заранее спасибо.

0

c++differential-equations runge-kutta

Решение

Другие решения

Других решений пока нет …

Источник

Accepted Answer

Обновить: Обнаружена одна ошибка: всегда компилировать с -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: выявлена вторая ошибка: в другом посте приведен код класса линейного вектора. Там, в отличие от прибавление (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])

как и ожидалось. вам понадобится отладчик или промежуточный вывод, чтобы найти, где вычисления идут неправильно.

0