Runge-Kutta duseno nihanje

Ko tudi učitelj ne more pomagati...
Post Reply
FilipGregor
Posts: 19
Joined: 1.12.2008 22:30

Runge-Kutta duseno nihanje

Post by FilipGregor » 4.8.2011 21:08

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})\)

User avatar
Aniviller
Posts: 7263
Joined: 15.11.2004 18:16

Re: Runge-Kutta duseno nihanje

Post by Aniviller » 5.8.2011 9:18

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
Posts: 19
Joined: 1.12.2008 22:30

Re: Runge-Kutta duseno nihanje

Post by FilipGregor » 5.8.2011 23:18

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.

User avatar
Aniviller
Posts: 7263
Joined: 15.11.2004 18:16

Re: Runge-Kutta duseno nihanje

Post by Aniviller » 6.8.2011 0:07

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
Posts: 19
Joined: 1.12.2008 22:30

Re: Runge-Kutta duseno nihanje

Post by FilipGregor » 6.8.2011 15:13

A tako. V praksi bi to tako zgledalo ne?

Code: Select all

// 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;
}

User avatar
Aniviller
Posts: 7263
Joined: 15.11.2004 18:16

Re: Runge-Kutta duseno nihanje

Post by Aniviller » 6.8.2011 16:27

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
Posts: 19
Joined: 1.12.2008 22:30

Re: Runge-Kutta duseno nihanje

Post by FilipGregor » 6.8.2011 16:32

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

sanej
Posts: 71
Joined: 25.8.2010 18:00

Re: Runge-Kutta duseno nihanje

Post by sanej » 13.6.2013 18:00

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

User avatar
Aniviller
Posts: 7263
Joined: 15.11.2004 18:16

Re: Runge-Kutta duseno nihanje

Post by Aniviller » 13.6.2013 23:33

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.

Post Reply