Runge-Kutta duseno nihanje

[FilipGregor]

Lepo pozdravljeni!

Z metodo Runge-Kutta 4. stopnje bi rad rešil diferencialno enačbo (v nadaljevanju d.e.) za dušeno nihanje. Vem da morem d.e. zapisati kot sistem dvem d.e. prvega reda. Problem imam s koeficientov k1, k2, k3 in k4. Prosil bi, za pomoč ker ne vem kako naj upoštevam tisti t+1/2 in y+(1/2)k1. Hvala!

d.e.

\(\ddot{x}=-\frac{b}{m}\dot{x}-\frac{k}{m}x\)

kot sistem d.e. prvega reda

\(\frac{dv}{dt} =-\frac{b}{m}\dot{x}-\frac{k}{m}x\)

\(\frac{dx}{dt}=v\)

Runge-Kutta metoda

\(\dot{y}=f(t,\ y)\)

\(y_{n+1}=y_{n}+\frac{1}{6}(k_{1}+2k_{2}+2k_{3}+k_{4})\)

\(k_{1}=h f(t_{n},\ y_{n})\)

\(k_{2}=h f(t_{n}+\frac{1}{2}h,\ y_{n}+\frac{1}{2}k_{1})\)

\(k_{3}=h f(t_{n}+\frac{1}{2}h,\ y_{n}+\frac{1}{2}k_{2})\)

\(k_{4}=h f(t_{n}+h,\ y_{n}+k_{3})\)

[Aniviller]

t+1/2 nima nobenega vpliva ker cas ne nastopa v tvoji diferencialni enacbi. Ostalo je pa takoj jasno, ce upostevas da je zdaj tvoj y vektor.
\(y=(x,v)\)
\(y'=(v,-b/m\cdot v-k/m\cdot x)\)
Zdaj imas zapisan ta y' kot funkcijo x in v (torej kot funkcijo y).
Zdaj bodo k1,k2,k3,... vektorji iz dveh komponent. Od tukaj naprej samo se izvedes postopek.

[FilipGregor]

Najlepša hvala za odgovor. Torej, če prav razumem se k-ji izračunajo tako?

\(k_{1}=\begin{bmatrix} v\\-\frac{b}{m}\ v-\frac{k}{m}\ x \end{bmatrix}\)

\(k_{2}=\begin{bmatrix} v+\frac{1}{2}\ v\\-\frac{b}{m}\ v-\frac{k}{m}\ x +\frac{1}{2}(-\frac{b}{m}\ v-\frac{k}{m}\ x) \end{bmatrix}\)

itd.

[Aniviller]

Ne, ti si naredil k2=f(y)+k1. k1 mora biti v argumentu. Torej
\(k_2=\begin{bmatrix}
v+h/2 (-b/m v-k/m x)\\
-b/m (v+h/2 (-b/m v-k/m x))-k/m(x+h/2 v)
\end{bmatrix}\)
Se pravi: kot "x" vstavis prvotni x plus x-komponento od k1 (krat tista polovicka koraka). In kot "v" vstavis prvotni v plus v-komponento od k1 (se pravi drugo).

Seveda se nikoli analiticno ne vstavlja, finta je v tem da to delas z numericnimi podatki. Poracunas stevilke za k1 in potem naprej.

[FilipGregor]

A tako. V praksi bi to tako zgledalo ne?

Koda: Izberi vse

// funkcija za izracun pospeska
double pos(double v, double x){
    return (-b/m*v - k/m*x);
}

// Runge-Kutta metoda
void RK4(double dt, double v, double x){
    k1v = pos(v,x)*dt;
    k1x = v * dt;
    k2v = pos(v + k1v/2 , x + k1x/2) * dt;
    k2x = (v + k1v/2)*dt;
    k3v = pos(v + k2v/2 , x + k2x/2) * dt;
    k3x = (v + k2v/2) * dt;
    k4v = pos(v + k3v , x + k3x) * dt;
    k4x = (v + k3x) * dt;

    vnp1 = v + (k1v + 2*k2v + 2*k3v + k4v) / 6.0;
    xnp1 = x + (k1x + 2*k2x + 2*k3x + k4x) / 6.0;
}

[Aniviller]

Tako nekako ja. Edino pri k4x ti nastopa notri k3x.

V splosnem, ce izraz za odvod x-a ne bi bil tako zelo enostaven, bi imel pa poleg "pos" se eno funkcijo za odvod komponente x. Tako da bi bili izrazi simetricni. Ponavadi numericne kniznjice, ki implementirajo te integratorje (v bistvu katerega koli), zahtevajo funkcijo, ki sprejme vektor in zapise rezultate v drug vektor. Recimo takole

void f(double t, double y[], double dy[]){
dy[0]=y[1];
dy[1]=-b/m*y[1]-k/m*y[0];
}
ce je y={x,v};

Tako da potem integratorja ne briga kako ti to racunas, pa tudi za mnozenje z dt sam poskrbi (ker ponavadi itak adaptivno spreminja korak glede na ocenjeno napako).

[FilipGregor]

OK, super! Najlepša hvala za pomoč!

[sanej]

Pozdravljeni!
Imam težave pri numeričnem računanju. Naloga pravi da je potrebno poiskati višino, hitrost in pospešek padalca.

\(\ddot{x} = \dot{x}^2 c \rho S/2m - g\) (1) to sem zapisal za višino

zanima me sledeče: reševal sem v matlabu in po gornji enačbi z metodo ode45 za reševanje de, sem višino izračunal pravilno, hitrost pa ki nastopa poleg pa je napačna. Navodilo pravi da je potrebno izračunati z Runge kutta reda 4. Ali je možno, da hitrost ne pride pravilno izračunana zaradi napake metode ?

Poleg tega pa še to: če izračunam hitrost iz (1), bi ta mogla biti identična hitrosti iz enačbe ali se motim ?? ker meni pridejo odstopanja velikosti 1m/s
\(\dot{v} = v^2 c \rho S/2 -g\)

Hvala za odgovore

[Aniviller]

Za takole enacbo bi morala biti kunge-kutta zelo robustna. Napake bi morale biti kvecjemu nekje na 7. decimalki ali kaj takega.

Sklepam, da resujes sistem
\(\dot{x}=v\)
\(\dot{v}=...\)
da dobis zraven poleg hitrosti se visino (ceprav v principu lahko integriras samo drugo enacbo in dobis potek hitrosti, to pa potem integriras se enkrat na katerikoli nacin, da dobis polozaj).

Pri drugi enacbi ti masa manjka, samo to je najbrz samo pozabljeno pri prepisovanju.