[C] Runge Kutta del quarto ordine
Buonasera, devo realizzare un programma in C che integri posizione e velocità di un oscillatore armonico con il metodo di Runge Kutta del quarto ordine. Io ho provato cosi con questa funzione:
Successivamente nel main ho:
A me sembra giusto questo codice ma ottengo risultati completamente sballati, per esempio per la posizone:
Qualcuno sa dirmi come mai?
double RungeKuttaFourth (double y, double dt, double dy) {
double y1, y2, y3, y4;
y1 = dy*dt;
y2 = (y + y1/2)*dt;
y3 = (y + y2/2)*dt;
y4 = (y + y3)*dt;
return y + 1.0/6.0*(y1 + 2*y2 + 2*y3 + y4);
}Successivamente nel main ho:
printf("\nHai scelto il metodo d'integrazione Runge-Kutta del quarto ordine");
for (i = 0; i <= N; i++) {
fprintf(f,"%lf %lf %lf\n", t, xANDv.x, xANDv.v);
xANDv.x = RungeKuttaFourth(xANDv.x,dt,xANDv.v);
xANDv.v = RungeKuttaFourth(xANDv.v,dt,-k/m*xANDv.x);
t+=dt;
}A me sembra giusto questo codice ma ottengo risultati completamente sballati, per esempio per la posizone:
Qualcuno sa dirmi come mai?
Risposte
Sia \(\displaystyle\frac{d\mathbf{y}}{dt} = \mathbf{f}( \mathbf{y}, t )\) un sistema di ODE che vuoi risovere. Allora RK4 si scrive come:
\(\displaystyle \mathbf{k}_1 = \mathbf{f}( \mathbf{y}_n, t_n ) \)
\(\displaystyle \mathbf{k}_2 = \mathbf{f}\biggl( \mathbf{y}_n + \frac{h\mathbf{k}_1}{2}, t_n + \frac{h}{2} \biggr) \)
\(\displaystyle \mathbf{k}_3 = \mathbf{f}\biggl( \mathbf{y}_n + \frac{h\mathbf{k}_2}{2}, t_n + \frac{h}{2} \biggr) \)
\(\displaystyle \mathbf{k}_4 = \mathbf{f}\biggl( \mathbf{y}_n + \mathbf{k}_3, t_n + h \biggr) \)
\(\displaystyle \mathbf{y}_{n+1} = \mathbf{y}_{n} + h\frac{\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4}{6} \)
Il tuo problema si scrive come:
\(\displaystyle\frac{d}{dt}\begin{pmatrix}x \\ v\end{pmatrix} = \mathbf{f}(x, v, t) = \begin{pmatrix}v \\ -\frac{k}{m}x\end{pmatrix}\)
A me non sembra che tu stia usando questa \(\mathbf{f}\).
\(\displaystyle \mathbf{k}_1 = \mathbf{f}( \mathbf{y}_n, t_n ) \)
\(\displaystyle \mathbf{k}_2 = \mathbf{f}\biggl( \mathbf{y}_n + \frac{h\mathbf{k}_1}{2}, t_n + \frac{h}{2} \biggr) \)
\(\displaystyle \mathbf{k}_3 = \mathbf{f}\biggl( \mathbf{y}_n + \frac{h\mathbf{k}_2}{2}, t_n + \frac{h}{2} \biggr) \)
\(\displaystyle \mathbf{k}_4 = \mathbf{f}\biggl( \mathbf{y}_n + \mathbf{k}_3, t_n + h \biggr) \)
\(\displaystyle \mathbf{y}_{n+1} = \mathbf{y}_{n} + h\frac{\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4}{6} \)
Il tuo problema si scrive come:
\(\displaystyle\frac{d}{dt}\begin{pmatrix}x \\ v\end{pmatrix} = \mathbf{f}(x, v, t) = \begin{pmatrix}v \\ -\frac{k}{m}x\end{pmatrix}\)
A me non sembra che tu stia usando questa \(\mathbf{f}\).
Forse è meglio se spiego cosa non va. Il problema è che \(x\) e \(v\) non sono indipendenti: non puoi chiamare rk4 due volte.
Scrivo \(\displaystyle \alpha = k/m \) e \(h = 2s\) per semplificare il testo. Noto inoltre che \(\displaystyle \mathbf{f} \) non dipende da \(t\) quindi:
\(\displaystyle \mathbf{f}( \mathbf{y} ) = \mathbf{f}\begin{pmatrix} x \\ v \end{pmatrix} = \begin{pmatrix} v \\ -\alpha x \end{pmatrix} \)
\(\displaystyle \mathbf{k}_1 = \mathbf{f}(\mathbf{y}) = \begin{pmatrix} v \\ -\alpha x \end{pmatrix} \)
\(\displaystyle \mathbf{k}_2 = \mathbf{f}(\mathbf{y} + s\mathbf{k}_1) = \mathbf{f}\begin{pmatrix} x + sv \\ v -\alpha s x \end{pmatrix} = \begin{pmatrix} v -\alpha s x \\ -\alpha (x + s v) \end{pmatrix} = \begin{pmatrix} v -\alpha s x \\ -\alpha x -\alpha s v \end{pmatrix} \)
\(\displaystyle \mathbf{k}_3 = \mathbf{f}(\mathbf{y} + s\mathbf{k}_2) = \mathbf{f}\begin{pmatrix} x + sv -\alpha s^2 x \\ v -\alpha s x -\alpha s^2 v \end{pmatrix} = \begin{pmatrix} v -\alpha s x -\alpha s^2 v \\ -\alpha ( x + sv -\alpha s^2 x) \end{pmatrix} = \begin{pmatrix} v -\alpha s x -\alpha s^2 v \\ -\alpha x -\alpha sv +\alpha^2 s^2 x \end{pmatrix} \)
e così anche per \(\mathbf{k}_4\).
Scrivo \(\displaystyle \alpha = k/m \) e \(h = 2s\) per semplificare il testo. Noto inoltre che \(\displaystyle \mathbf{f} \) non dipende da \(t\) quindi:
\(\displaystyle \mathbf{f}( \mathbf{y} ) = \mathbf{f}\begin{pmatrix} x \\ v \end{pmatrix} = \begin{pmatrix} v \\ -\alpha x \end{pmatrix} \)
\(\displaystyle \mathbf{k}_1 = \mathbf{f}(\mathbf{y}) = \begin{pmatrix} v \\ -\alpha x \end{pmatrix} \)
\(\displaystyle \mathbf{k}_2 = \mathbf{f}(\mathbf{y} + s\mathbf{k}_1) = \mathbf{f}\begin{pmatrix} x + sv \\ v -\alpha s x \end{pmatrix} = \begin{pmatrix} v -\alpha s x \\ -\alpha (x + s v) \end{pmatrix} = \begin{pmatrix} v -\alpha s x \\ -\alpha x -\alpha s v \end{pmatrix} \)
\(\displaystyle \mathbf{k}_3 = \mathbf{f}(\mathbf{y} + s\mathbf{k}_2) = \mathbf{f}\begin{pmatrix} x + sv -\alpha s^2 x \\ v -\alpha s x -\alpha s^2 v \end{pmatrix} = \begin{pmatrix} v -\alpha s x -\alpha s^2 v \\ -\alpha ( x + sv -\alpha s^2 x) \end{pmatrix} = \begin{pmatrix} v -\alpha s x -\alpha s^2 v \\ -\alpha x -\alpha sv +\alpha^2 s^2 x \end{pmatrix} \)
e così anche per \(\mathbf{k}_4\).
Forse ho capito cosa intendi, una volta calcolato y1 sia per posizione che velocità per calcolare y2 devo usare i nuovi x e v ricavati da y1?
Ci sono due problemi nel tuo codice:
1. Come già spiegato da @vict85, lo stato del sistema comprende sia la posizione che la velocità e questi valori devono essere entrambi passati alla funzione.
2. La derivata viene calcolata quattro volte nel metodo e non è possibile passare semplicemente il suo valore al tempo t alla funzione.
Ho scritto un piccolo codice di esempio per mostrare come implementare questo metodo (per quanto in pratica il codice sia di solito più complesso e generale di quello che ho scritto). Visualizzando i valori ottenuti mi ritrovo effettivamente una sinusoide come ci si aspetterebbe in questo caso.
Nota come tutto lo stato sia passato alle funzioni e non solo posizione o velocità separatamente e che la funzione che calcola la derivata sia chiamata più volte ad ogni iterazione.
1. Come già spiegato da @vict85, lo stato del sistema comprende sia la posizione che la velocità e questi valori devono essere entrambi passati alla funzione.
2. La derivata viene calcolata quattro volte nel metodo e non è possibile passare semplicemente il suo valore al tempo t alla funzione.
Ho scritto un piccolo codice di esempio per mostrare come implementare questo metodo (per quanto in pratica il codice sia di solito più complesso e generale di quello che ho scritto). Visualizzando i valori ottenuti mi ritrovo effettivamente una sinusoide come ci si aspetterebbe in questo caso.
#include <math.h>
#include <stdio.h>
// Dynamic system parameters
const double k = 2.0;
const double m = 3.0;
// Dynamic system state (position and velocity)
typedef struct state_s {
double x;
double v;
} state_s;
// Implements a fused multiply-add operation on the state.
// returns `s0 + (t * s1)`
state_s state_fma(state_s s0, double t, state_s s1)
{
return (state_s){
.x = fma(t, s1.x, s0.x),
.v = fma(t, s1.v, s0.v)
};
}
// Derivative of the dynamic system
state_s f(double t, state_s state)
{
return (state_s){
.x = state.v,
.v = - (k / m) * state.x
};
}
// Runge-Kutta of 4th order
state_s rk4(state_s state, double t, double dt)
{
double h = 0.5 * dt;
state_s k1 = f(t, state);
state_s k2 = f(t + h, state_fma(state, h, k1));
state_s k3 = f(t + h, state_fma(state, h, k2));
state_s k4 = f(t + dt, state_fma(state, dt, k3));
double a = dt / 6.0;
state = state_fma(state, a, k1);
state = state_fma(state, 2.0 * a, k2);
state = state_fma(state, 2.0 * a, k3);
state = state_fma(state, a, k4);
return state;
}
// Create table of values
int main()
{
FILE *fout = fopen("sim.csv", "w");
state_s state = (state_s){ .x = 0.0, .v = 1.0 };
double dt = 0.1;
double t = 0.0;
for (int i = 0; i < 100; ++i) {
fprintf(fout, "%f %f\n", state.x, state.v);
state = rk4(state, t, dt);
t += dt;
}
fprintf(fout, "%f %f\n", state.x, state.v);
fclose(fout);
}
Nota come tutto lo stato sia passato alle funzioni e non solo posizione o velocità separatamente e che la funzione che calcola la derivata sia chiamata più volte ad ogni iterazione.
Non mi è chiarissimo quello che hai fatto. Ho pensato di rifare la funzione cosi, pensando all'indipendenza di x e v:
La situazione migliora un pò però ancora non ci siamo
xv OscillatorRK4 (xv XV_old, double *t, double *dE, double dt, double m, double k, double E0) {
xv XV_new;
xv xv1, xv2, xv3, xv4;
xv1.x = XV_old.v*dt;
xv1.v = -k/m*XV_old.x*dt;
xv2.x = (XV_old.v + xv1.v/2.0)*dt;
xv2.v = -k/m*(XV_old.x + xv1.x/2.0)*dt;
xv3.x = (XV_old.v + xv2.v/2.0)*dt;
xv3.v = -k/m*(XV_old.x + xv2.x/2.0)*dt;
xv4.x = (XV_old.v + xv3.v)*dt;
xv4.v = -k/m*(XV_old.x + xv3.x/2.0)*dt;
XV_new.x = XV_old.x + (xv1.x + 2*xv2.x + 2* xv3.x + 2*xv4.x)/6.0;
XV_new.v = XV_old.v + (xv1.v + 2*xv2.v + 2* xv3.v + 2*xv4.v)/6.0;
*t += dt;
*dE = (Energy(m,k,XV_new.v,XV_new.x) - E0)/E0;
return XV_new;
}La situazione migliora un pò però ancora non ci siamo
Come ti abbiamo già spiegato, il metodo Runge-Kutta del quarto ordine permette di integrare un sistema alle derivate parziali ordinario del primo ordine nella forma:
\[ \frac{dx}{dt} = \dot{x} = f(t, x) \]
con una condizione iniziale del tipo \(x(t_0) = x_0\). È un metodo ricorsivo che calcola un nuovo valore \(x_{n+1}\) a partire dal valore \(x_n\) e dalla funzione \(f(t, x)\). Per prima cosa ti devi calcolare
\[ \begin{align*}
k_1 &= f(t_n, x_n) \\
k_2 &= f(t_n + \Delta t/2, x_n + (\Delta t/2) \, k_1) \\
k_3 &= f(t_n + \Delta t/2, x_n + (\Delta t/2) \, k_2) \\
k_4 &= f(t_n + \Delta t, x_n + \Delta t \, k_3)
\end{align*} \]
e usa questi valori per calcolare infine
\[ x_{n+1} = x_n + \frac{h\,(k_1 + 2\,k_2 + 2\,k_3 + k_4)}{6} \]
Devi valutare il valore di \(f\) quattro volte nel metodo. Il codice che ti avevo scritto calcolava infatti il valore di \(f\) quattro volte per trovare i vari \(k_i\) per poi sommare tra di loro i risultati. Che cosa non ti era chiaro esattamente del codice?
È poi invece possibile scrivere un codice il cui unico obiettivo sia quello di risolvere il tuo problema specifico ma in effetti perde completamente la forma che ho descritto sopra perché diventa una semplice moltiplicazione matriciale come ho scritto più in basso. Nel seguito ti faccio vedere come scrivere una versione ottimale del metodo per il tuo problema specifico.
Nel tuo caso vuoi risolvere la seguente equazione differenziale del secondo ordine:
\[ \frac{d^2 x}{dt^2} = \ddot{x} = f(t, x) = - \alpha\,x. \]
Il metodo classico per risolvere tale equazione è quella di riscrivere l'equazione come un sistema in cui le derivate intermedie sono sostituite da variabili intermedie per cui dovrai risolvere
\[ \begin{cases}
dx/dt = \dot{x} = v \\
dv/dt = \dot{v} = f(t, x) = - \alpha\,x
\end{cases} \]
Si tratta di un sistema differenziale lineare nella forma:
\[ \dot{Z} =
\begin{bmatrix}
\dot{x} \\
\dot{v}
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 \\
- \alpha & 0
\end{bmatrix}
\begin{bmatrix}
x \\
v
\end{bmatrix}
=
J\,Z
\]
In questo caso hai che
\[
\begin{align*}
K_1 &= J\,Z_n \\
K_2 &= J\,(Z_n + (\Delta t/2)\,K_1) = \bigl(J + (\Delta t / 2)\,J^2\bigr)\,Z_n \\
K_3 &= J\,(Z_n + (\Delta t/2)\,K_2) = \bigl(J + (\Delta t / 2)\,J^2 + (\Delta t / 2)^2\,J^3\bigr)\,Z_n \\
K_4 &= J\,(Z_n + \Delta t\;K_3) = \bigl(J + \Delta t\,J^2 + (\Delta t^2 / 2)\,J^3 + (\Delta t^3 / 4)\,J^4 \bigr)\,Z_n \\
\end{align*}
\]
A questo punto puoi scrivere (se non ho sbagliato i calcoli)
\[ Z_{n+1} = \left(I + \Delta t\,J + \frac{\Delta t^2\,J^2}{2} + \frac{\Delta t^3\,J^3}{6} + \frac{\Delta t^4\,J^4}{24} \right)\, Z_n \]
Siccome tutta quella espressione è costante puoi calcolarla in anticipo ottenendo (sempre se non ho sbagliato i calcoli)
\[
\begin{bmatrix}
\alpha^2 \, \Delta t^4 \, / 24 - \alpha \, \Delta t^2 \, / 2 + 1 & - \alpha\,\Delta t^3 / 6 + \Delta t \\
\alpha^2 \, \Delta t^3 / 6 - \alpha\,\Delta t & \alpha^2\,\Delta t^4\,/24 - \alpha\,\Delta t^2 / 2 + 1
\end{bmatrix}
\]
Il seguente codice fa uso di questa ottimizzazione:
Nota come ad ogni iterazione il metodo si riduca a quattro moltiplicazioni e due somme..
\[ \frac{dx}{dt} = \dot{x} = f(t, x) \]
con una condizione iniziale del tipo \(x(t_0) = x_0\). È un metodo ricorsivo che calcola un nuovo valore \(x_{n+1}\) a partire dal valore \(x_n\) e dalla funzione \(f(t, x)\). Per prima cosa ti devi calcolare
\[ \begin{align*}
k_1 &= f(t_n, x_n) \\
k_2 &= f(t_n + \Delta t/2, x_n + (\Delta t/2) \, k_1) \\
k_3 &= f(t_n + \Delta t/2, x_n + (\Delta t/2) \, k_2) \\
k_4 &= f(t_n + \Delta t, x_n + \Delta t \, k_3)
\end{align*} \]
e usa questi valori per calcolare infine
\[ x_{n+1} = x_n + \frac{h\,(k_1 + 2\,k_2 + 2\,k_3 + k_4)}{6} \]
Devi valutare il valore di \(f\) quattro volte nel metodo. Il codice che ti avevo scritto calcolava infatti il valore di \(f\) quattro volte per trovare i vari \(k_i\) per poi sommare tra di loro i risultati. Che cosa non ti era chiaro esattamente del codice?
È poi invece possibile scrivere un codice il cui unico obiettivo sia quello di risolvere il tuo problema specifico ma in effetti perde completamente la forma che ho descritto sopra perché diventa una semplice moltiplicazione matriciale come ho scritto più in basso. Nel seguito ti faccio vedere come scrivere una versione ottimale del metodo per il tuo problema specifico.
Nel tuo caso vuoi risolvere la seguente equazione differenziale del secondo ordine:
\[ \frac{d^2 x}{dt^2} = \ddot{x} = f(t, x) = - \alpha\,x. \]
Il metodo classico per risolvere tale equazione è quella di riscrivere l'equazione come un sistema in cui le derivate intermedie sono sostituite da variabili intermedie per cui dovrai risolvere
\[ \begin{cases}
dx/dt = \dot{x} = v \\
dv/dt = \dot{v} = f(t, x) = - \alpha\,x
\end{cases} \]
Si tratta di un sistema differenziale lineare nella forma:
\[ \dot{Z} =
\begin{bmatrix}
\dot{x} \\
\dot{v}
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 \\
- \alpha & 0
\end{bmatrix}
\begin{bmatrix}
x \\
v
\end{bmatrix}
=
J\,Z
\]
In questo caso hai che
\[
\begin{align*}
K_1 &= J\,Z_n \\
K_2 &= J\,(Z_n + (\Delta t/2)\,K_1) = \bigl(J + (\Delta t / 2)\,J^2\bigr)\,Z_n \\
K_3 &= J\,(Z_n + (\Delta t/2)\,K_2) = \bigl(J + (\Delta t / 2)\,J^2 + (\Delta t / 2)^2\,J^3\bigr)\,Z_n \\
K_4 &= J\,(Z_n + \Delta t\;K_3) = \bigl(J + \Delta t\,J^2 + (\Delta t^2 / 2)\,J^3 + (\Delta t^3 / 4)\,J^4 \bigr)\,Z_n \\
\end{align*}
\]
A questo punto puoi scrivere (se non ho sbagliato i calcoli)
\[ Z_{n+1} = \left(I + \Delta t\,J + \frac{\Delta t^2\,J^2}{2} + \frac{\Delta t^3\,J^3}{6} + \frac{\Delta t^4\,J^4}{24} \right)\, Z_n \]
Siccome tutta quella espressione è costante puoi calcolarla in anticipo ottenendo (sempre se non ho sbagliato i calcoli)
\[
\begin{bmatrix}
\alpha^2 \, \Delta t^4 \, / 24 - \alpha \, \Delta t^2 \, / 2 + 1 & - \alpha\,\Delta t^3 / 6 + \Delta t \\
\alpha^2 \, \Delta t^3 / 6 - \alpha\,\Delta t & \alpha^2\,\Delta t^4\,/24 - \alpha\,\Delta t^2 / 2 + 1
\end{bmatrix}
\]
Il seguente codice fa uso di questa ottimizzazione:
#include <math.h>
#include <stdio.h>
// Dynamic system parameters
const double k = 2.0;
const double m = 3.0;
const double alpha = k / m;
const double dt = 0.1;
const double tmax = 10.0;
const double x0 = 5.0;
// Runge-Kutta constants
const double alpha2 = alpha * alpha;
const double dt2 = dt * dt;
const double dt3 = dt2 * dt;
const double dt4 = dt2 * dt2;
const double M00 = alpha2 * dt4 / 24.0 - alpha * dt2 / 2.0 + 1.0;
const double M01 = - alpha * dt3 / 6.0 + dt;
const double M10 = alpha2 * dt3 / 6.0 - alpha * dt;
const double M11 = alpha2 * dt4 / 24 - alpha * dt2 / 2 + 1.0;
// Runge-Kutta of 4th order
void rk4(double old_state[2], double dt, double new_state[2])
{
new_state[0] = M00 * old_state[0] + M01 * old_state[1];
new_state[1] = M10 * old_state[0] + M11 * old_state[1];
}
// Create table of values
int main()
{
FILE *fout = fopen("sim.csv", "w");
double t = 0.0;
double Z[2] = {x0, 0.0};
fprintf(fout, "%f %f %f %f\n", t, Z[0], Z[1], 0.0);
while (t < (tmax - 0.5 * dt)) {
double newZ[2] = {0.0, 0.0};
const double freq = sqrtf(alpha);
double sol = x0 * cos(freq * (t + dt));
rk4(Z, dt, newZ);
t += dt;
Z[0] = newZ[0];
Z[1] = newZ[1];
fprintf(fout, "%f %f %f %f\n", t, Z[0], Z[1], fabs(Z[0] - sol));
}
fclose(fout);
}
Nota come ad ogni iterazione il metodo si riduca a quattro moltiplicazioni e due somme..
Grazie mille ho capito ora. Non perchè il tuo non sia corretto, anzi è molto più pulito del mio, ma preferisco produrre qualcosa per conto mio 
E nel main:
Questo codice fa anche altre cose, come ad esempio vedere come varia l'energia.
Per la posizione ottengo questi risultati:
Mi sembra sia accettabile, voi che dite?
xv OscillatorRK4 (xv XV_old, double *t, double *dE, double dt, double m, double k, double E0) {
xv XV_new;
xv xv1, xv2, xv3, xv4;
xv1.x = XV_old.v*dt;
xv1.v = -k/m*XV_old.x*dt;
xv2.x = (XV_old.v + xv1.v/2)*dt;
xv2.v = -k/m*(XV_old.x+xv1.x/2)*dt;
xv3.x = (XV_old.v + xv2.v/2)*dt;
xv3.v = -k/m*(XV_old.x+xv2.x/2)*dt;
xv4.x = (XV_old.v + xv3.v)*dt;
xv4.v = -k/m*(XV_old.x+xv3.x)*dt;
XV_new.x = XV_old.x + (xv1.x + 2*xv2.x + 2*xv3.x + xv4.x)/6;
XV_new.v = XV_old.v + (xv1.v + 2*xv2.v + 2*xv3.v + xv4.v)/6;
*t += dt;
*dE = (Energy(m,k,XV_new.v,XV_new.x) - E0)/E0;
return XV_new;
}E nel main:
else if (metodo == RUNGEKUTTA4) {
printf("\nHai scelto il metodo d'integrazione Runge-Kutta del quarto ordine");
for (i = 0; i <= N; i++) {
fprintf(f,"%lf %lf %lf %lf\n", t, xANDv.x, xANDv.v, dE);
xANDv = OscillatorRK4(xANDv,&t,&dE,dt,m,k,E0);
}
}Questo codice fa anche altre cose, come ad esempio vedere come varia l'energia.
Per la posizione ottengo questi risultati:
Mi sembra sia accettabile, voi che dite?
È un bene implementare per conto proprio, ma si deve anche cercare di imparare dai codici altrui. Per esempio, l'ultimo messaggio di apatriarca ti mostra come ottimizzare un algoritmo generico ad un problema specifico. Anche se non ha senso impegnarsi tanto per ogni algoritmo. Il primo per esempio mostra l'uso della funzione fma (fused multiplied and add).
Nota che il tuo codice non è troppo leggibile. Sinceramente userei un nome più comprensibile di [inline]xv[/inline] per il tipo e, in C, tendo a non usare il typedef per le strutture (ma questo è un preferenza personale). Se non preservi lo struct all'inizio, una scelta classica è [inline]stato_s[/inline]. Inoltre non chiamerei tutte le variabile di tipo [inline]xv[/inline] come [inline]xv_qualcosa[/inline]. Per esempio, per le varie [inline]Kn[/inline] userei semplicemente [inline]kn[/inline].
Invece di fare [inline]XV_old.x[/inline] e [inline]XV_old.v[/inline] ogni volta, puoi inoltre creare le variabili locali [inline]x[/inline] e [inline]v[/inline]. Inoltre precalcolerei [inline]-k/m[/inline].
Che editor usi? Dovresti pensare di usare un qualche programma per formattare il codice. Io uso clang-format che può essere usato da quasi tutti i principali editor. Quelli che non lo hanno, hanno qualcos'altro.
Nota che il tuo codice non è troppo leggibile. Sinceramente userei un nome più comprensibile di [inline]xv[/inline] per il tipo e, in C, tendo a non usare il typedef per le strutture (ma questo è un preferenza personale). Se non preservi lo struct all'inizio, una scelta classica è [inline]stato_s[/inline]. Inoltre non chiamerei tutte le variabile di tipo [inline]xv[/inline] come [inline]xv_qualcosa[/inline]. Per esempio, per le varie [inline]Kn[/inline] userei semplicemente [inline]kn[/inline].
Invece di fare [inline]XV_old.x[/inline] e [inline]XV_old.v[/inline] ogni volta, puoi inoltre creare le variabili locali [inline]x[/inline] e [inline]v[/inline]. Inoltre precalcolerei [inline]-k/m[/inline].
Che editor usi? Dovresti pensare di usare un qualche programma per formattare il codice. Io uso clang-format che può essere usato da quasi tutti i principali editor. Quelli che non lo hanno, hanno qualcos'altro.
È giusto che tu scriva il codice di conto tuo. I miei esempi avevano solo lo scopo di mostrare in codice quello che volevo trasmettere nel testo del messaggio. Siccome esiste un metodo semplice per ottenere la soluzione analitica, ti consiglio di calcolarti anche l'errore commesso. Perché sia corretto dovrebbe crescere secondo un determinato limite asintotico.
Una delle ragioni per cui ti ho scritto i due codici era mostrare quelli che credo siano i due principali modi per scrivere una funzione del genere:
1. Come una funzione generica che può essere usata per diversi problemi cambiando la funzione o la dimensione del problema come nel primo esempio che ti ho scritto. Volessi cambiare il problema con qualcosa di completamente diverso potrei cambiare le funzioni di contorno, ma la funzione del metodo sarebbe sempre uguale.
2. Ottimizzata per il problema corrente e in questo caso può valer la pena di ottimizzare al massimo come ho fatto nel secondo caso.
Il tuo codice è un po' una via di mezzo. Non è infatti immediato modificarlo se si cambia il problema anche di poco (considera per esempio se aggiungi uno smorzatore*—come cambierebbe il codice?). Ma è tuttavia anche molto lontano dall'essere una versione ottimizzata per il problema specifico. Il mio consiglio è quello di aggiungere qualche funzione di supporto e di aggiungere qualche spazio ulteriore.
* Con aggiungere uno smorzatore intendo dire risolvere l'equazione differenziale
\[\ddot{x} + 2\,\zeta\,\omega\,\dot{x} + \omega^2\,x = 0 \]
dove \(\omega = \sqrt{\frac{k}{m}}\), con \(k\) la costante della molla e \(m\) la massa dell´oggetto, e \(\zeta = \frac{c}{2\,\sqrt{m\,k}}\) con \(c\) la costante di smorzamento. Hai aggiunto insomma una forza uguale a \(- c\,\dot{x}\) al sistema. Anche in questo caso esiste una soluzione analitica al problema anche se è un po' più complicata.
Una delle ragioni per cui ti ho scritto i due codici era mostrare quelli che credo siano i due principali modi per scrivere una funzione del genere:
1. Come una funzione generica che può essere usata per diversi problemi cambiando la funzione o la dimensione del problema come nel primo esempio che ti ho scritto. Volessi cambiare il problema con qualcosa di completamente diverso potrei cambiare le funzioni di contorno, ma la funzione del metodo sarebbe sempre uguale.
2. Ottimizzata per il problema corrente e in questo caso può valer la pena di ottimizzare al massimo come ho fatto nel secondo caso.
Il tuo codice è un po' una via di mezzo. Non è infatti immediato modificarlo se si cambia il problema anche di poco (considera per esempio se aggiungi uno smorzatore*—come cambierebbe il codice?). Ma è tuttavia anche molto lontano dall'essere una versione ottimizzata per il problema specifico. Il mio consiglio è quello di aggiungere qualche funzione di supporto e di aggiungere qualche spazio ulteriore.
* Con aggiungere uno smorzatore intendo dire risolvere l'equazione differenziale
\[\ddot{x} + 2\,\zeta\,\omega\,\dot{x} + \omega^2\,x = 0 \]
dove \(\omega = \sqrt{\frac{k}{m}}\), con \(k\) la costante della molla e \(m\) la massa dell´oggetto, e \(\zeta = \frac{c}{2\,\sqrt{m\,k}}\) con \(c\) la costante di smorzamento. Hai aggiunto insomma una forza uguale a \(- c\,\dot{x}\) al sistema. Anche in questo caso esiste una soluzione analitica al problema anche se è un po' più complicata.
Si capisco cosa intendi, purtroppo non sono ancora in grado di realizzare programmi "portatili" e ho sempre da modificare qualche cosa ad ogni piccola variazione. Comunque ho provato a stimare l'errore sull'energia al variare di dt e mi viene una curva del tipo $f(x) = ax^5 $ o giù di li. Il prof mi ricordo disse che per l'oscillatore e per piccoli dt era un metodo dell'ordine 5 quindi credo ci posso stare. L'oscillatore smorzato e con forzante l'ho provato già con RK II. però come dici tu ho dovuto passare quei 10 minuti a sistemare delle cosette. Comunque grazie dei consigli, proverò a migliorare nei prossimi programmi.
Portabilità si riferisce al trasferire il codice su un ambiente di esercuzione diverso (per esempio un diverso sistema operativo). Penso che tu intedessi che non sai realizzare programmi sufficientemente generici da poter essere riutilizzato per problemi diversi. Comunque anche quella è una questione di esercizio.
Si hai ragione, ho fatto confusione.
Aggiungerei solamente che con RK4, nel caso del pendolo semplice (cioè senza forzante), non avrai conservazione dell'energia anche a livello della soluzione discreta. Ossia, se plotti in un grafico il tempo e in ordinata l'energia del sistema $E(t) = K+V = \frac{1}{2}m v(t)^2 + V(t) $ , osserverai un "gradino"!
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo