Soluzione approssimata, equazione differenziale runge-kutta
Ero indeciso se creare il post su analisi matematica o qui, mi scuso se ho sbagliato.
Ho un problema nella comprensione dell'algoritmo di runge-kutta.
Devo approssimare la soluzione di questo sistema di equazioni differenziali attraverso runge-kutta:
[tex]\left\{
\begin{array}{lr}
\frac{d\theta}{dt} = \omega = f_1(t, \omega) \\
\frac{d\omega}{dt} = A * sin(\theta) + B*v_a - B * \omega = f_2(t, \omega, \theta)
\end{array}}[/tex]
L'algoritmo che sto studiando è RK4:
Per risolvere [tex]\frac{dy}{dt} = f(t, x)[/tex]
[tex]h = \frac{x_{finale} - x_{iniziale}}{N}[/tex] è il passo di integrazione
[tex]k_1 = f(x_n, y_n)[/tex]
[tex]k_2 = f(x_n + \frac{h}{2}, y_n +\frac{k_1*h}{2})[/tex]
[tex]k_3 = f(x_n + \frac{h}{2}, y_n + \frac{k_2*h}{2})[/tex]
[tex]k_4 = f(x_n + h, y_n + k_3*h)[/tex]
[tex]y_{n+1} = y_n + \frac{k_1}{6} + \frac{k_2}{3} + \frac{k_3}{3} + \frac{k_4}{6} + O(h^5)[/tex]
Ma non riesco ad applicarlo perchè $f2$ ha più variabili.
Potreste darmi una mano?
Grazie
Ho un problema nella comprensione dell'algoritmo di runge-kutta.
Devo approssimare la soluzione di questo sistema di equazioni differenziali attraverso runge-kutta:
[tex]\left\{
\begin{array}{lr}
\frac{d\theta}{dt} = \omega = f_1(t, \omega) \\
\frac{d\omega}{dt} = A * sin(\theta) + B*v_a - B * \omega = f_2(t, \omega, \theta)
\end{array}}[/tex]
L'algoritmo che sto studiando è RK4:
Per risolvere [tex]\frac{dy}{dt} = f(t, x)[/tex]
[tex]h = \frac{x_{finale} - x_{iniziale}}{N}[/tex] è il passo di integrazione
[tex]k_1 = f(x_n, y_n)[/tex]
[tex]k_2 = f(x_n + \frac{h}{2}, y_n +\frac{k_1*h}{2})[/tex]
[tex]k_3 = f(x_n + \frac{h}{2}, y_n + \frac{k_2*h}{2})[/tex]
[tex]k_4 = f(x_n + h, y_n + k_3*h)[/tex]
[tex]y_{n+1} = y_n + \frac{k_1}{6} + \frac{k_2}{3} + \frac{k_3}{3} + \frac{k_4}{6} + O(h^5)[/tex]
Ma non riesco ad applicarlo perchè $f2$ ha più variabili.
Potreste darmi una mano?
Grazie
Risposte
Stai lavorando con un sistema di equazioni differenziali del primo ordine, dunque nella forma
[tex]$\begin{cases} \underline y'(t) = \underline F(t, \underline y(t)) \\ \underline y(t_0) = \underline y_0\end{cases}$[/tex]
Nel tuo caso [tex]$\underline y(t) = \begin{bmatrix}\theta(t) \\ \omega(t)\end{bmatrix}$[/tex] e [tex]$F(t,\underline y(t)) = \begin{bmatrix}\omega(t) \\ A*\sin(\theta(t)) + B*v_a - B*\omega(t)\end{bmatrix}$[/tex]
Dunque rispetto al problema generale avrai
[tex]$\underline k_1=\underline F(x_n, \underline y_n)$[/tex]
[tex]$\underline k_2=\underline F(x_n + \frac{h}{2}, \underline y_n + \underline k_1 \frac{h}{2})$[/tex]
[tex]$\underline k_3=\underline F(x_n + \frac{h}{2}, \underline y_n + \underline k_2 \frac{h}{2})$[/tex]
[tex]$\underline k_4=\underline F(x_n + h, \underline y_n + \underline k_3 h)$[/tex]
[tex]$y_{n+1} = y_n + \frac{h}{6}(\underline k_1 + 2 \underline k_2 + 2 \underline k_3 + \underline k_4 ) + O(h^5)$[/tex]
[tex]$\begin{cases} \underline y'(t) = \underline F(t, \underline y(t)) \\ \underline y(t_0) = \underline y_0\end{cases}$[/tex]
Nel tuo caso [tex]$\underline y(t) = \begin{bmatrix}\theta(t) \\ \omega(t)\end{bmatrix}$[/tex] e [tex]$F(t,\underline y(t)) = \begin{bmatrix}\omega(t) \\ A*\sin(\theta(t)) + B*v_a - B*\omega(t)\end{bmatrix}$[/tex]
Dunque rispetto al problema generale avrai
[tex]$\underline k_1=\underline F(x_n, \underline y_n)$[/tex]
[tex]$\underline k_2=\underline F(x_n + \frac{h}{2}, \underline y_n + \underline k_1 \frac{h}{2})$[/tex]
[tex]$\underline k_3=\underline F(x_n + \frac{h}{2}, \underline y_n + \underline k_2 \frac{h}{2})$[/tex]
[tex]$\underline k_4=\underline F(x_n + h, \underline y_n + \underline k_3 h)$[/tex]
[tex]$y_{n+1} = y_n + \frac{h}{6}(\underline k_1 + 2 \underline k_2 + 2 \underline k_3 + \underline k_4 ) + O(h^5)$[/tex]
Grazie per la risposta.
Se ho capito bene devo iterare il processo per ogni funzione (quindi calcolare 2 $k_i$ e 2 $y_{i+1}$ ad ogni passo.
Ho provato a scrivere il codice in matlab:
Questa è la funzione che applica l'algoritmo:
Queste due sono le funzioni del sistema:
Che chiamo in questo modo (il penultimo parametro non serve e devo eliminarlo):
Ottengo un andamento oscillante a cui non so dare interpretazione ma non capisco se sbaglio l'intervallo, se sbaglio l'algoritmo o altro
Provo a ripensare a tutto quello che ho scritto.
Se ho scritto qualcosa di sbagliato, di cui non mi sono accorto, potreste avvertirmi?
Se ho capito bene devo iterare il processo per ogni funzione (quindi calcolare 2 $k_i$ e 2 $y_{i+1}$ ad ogni passo.
Ho provato a scrivere il codice in matlab:
Questa è la funzione che applica l'algoritmo:
Queste due sono le funzioni del sistema:
Che chiamo in questo modo (il penultimo parametro non serve e devo eliminarlo):
myRunge('mffunc', 'g', 0, 8, 500, 0, 0)
Ottengo un andamento oscillante a cui non so dare interpretazione ma non capisco se sbaglio l'intervallo, se sbaglio l'algoritmo o altro
Provo a ripensare a tutto quello che ho scritto.
Se ho scritto qualcosa di sbagliato, di cui non mi sono accorto, potreste avvertirmi?
No no, puoi fare tutto manipolando direttamente i vettori
Questa function realizza il metodo RK4, il formato è quello più simile alle function presenti in matlab:
- odefun rappresenta la funzione vettoriale nelle variabili [tex]$t$[/tex] e [tex]$\underline y$[/tex]
- tspan rappresenta gli estremi della variabile indipendente
- y0 è il vettore che rappresenta le condizioni iniziali
- Nh è il numero di intervallini in cui si suddivide l'intervallo considerato della variaible temporale
Si costruisce una matrice $u$ di dimensione $n \times T$, dove $n$ è il numero di componenti della variabile $y$ e $T$ è il numero degli istanti temporali considerati.
Si inizializza la prima colonna della matrice con $y_0$ condizione inziale, a questo punto si itera sui rimanenti istanti temporali, e si applica il metodo RK considerando i vettori nei vari istanti temporali.
Alla fine si traspone $u$ per mantenere la convenzione che il risultato contiene nelle colonne le soluzioni di ciascuna componente.
La tua funzione $F$ da passare dovrebbe essere ad esempio
function [t, u] = rk4(odefun, tspan, y0, Nh, varargin)
Nh=ceil(Nh);
t=linspace(tspan(1),tspan(2), Nh + 1);
u = zeros(length(y0), length(t));
u(:,1) = y0;
h=t(2)-t(1);
for i=2:length(t)
K1 = odefun(t(i-1), u(:,i-1));
K2 = odefun(t(i-1) + h/2, u(:,i-1) + h/2*K1);
K3 = odefun(t(i-1) + h/2, u(:,i-1) + h/2*K2);
K4 = odefun(t(i), u(:,i-1) + h*K3);
u(:,i) = u(:,i-1) + h/6*(K1 + 2*K2 + 2*K3 + K4);
end;
u=u';
Questa function realizza il metodo RK4, il formato è quello più simile alle function presenti in matlab:
- odefun rappresenta la funzione vettoriale nelle variabili [tex]$t$[/tex] e [tex]$\underline y$[/tex]
- tspan rappresenta gli estremi della variabile indipendente
- y0 è il vettore che rappresenta le condizioni iniziali
- Nh è il numero di intervallini in cui si suddivide l'intervallo considerato della variaible temporale
Si costruisce una matrice $u$ di dimensione $n \times T$, dove $n$ è il numero di componenti della variabile $y$ e $T$ è il numero degli istanti temporali considerati.
Si inizializza la prima colonna della matrice con $y_0$ condizione inziale, a questo punto si itera sui rimanenti istanti temporali, e si applica il metodo RK considerando i vettori nei vari istanti temporali.
Alla fine si traspone $u$ per mantenere la convenzione che il risultato contiene nelle colonne le soluzioni di ciascuna componente.
La tua funzione $F$ da passare dovrebbe essere ad esempio
F=@(t,y)[y(2); A*sin(y(1))+B*v_a-B*y(2)];
Grazie, provo ad utilizzare il codice che hai postato applicandola al mio caso.
La funzione F di cui parli è quella che passerò come odefun giusto?
La funzione F di cui parli è quella che passerò come odefun giusto?
Ho provato a risolvere l'equazione come ti ho suggerito, e il risultato è un andamento oscillate.
http://skadotnet.altervista.org/test.m
http://skadotnet.altervista.org/rk4.m
http://skadotnet.altervista.org/test.m
http://skadotnet.altervista.org/rk4.m
si anche io ho avuto lo stesso risultato.
Però se aumenti il lasso di tempo la funzione si assesta (dovrebbe essere un punto di equilibrio)
Proverò con altri valori in base allo studio del sistema che ho fatto.
Grazie mille. Sono commosso, sei il primo che mi risponde.
Però se aumenti il lasso di tempo la funzione si assesta (dovrebbe essere un punto di equilibrio)
Proverò con altri valori in base allo studio del sistema che ho fatto.
Grazie mille. Sono commosso, sei il primo che mi risponde.
ciao a tutti...scusate se mi intrometto...
io devo effettuare l'integrazione con runge-kutta
dO/dt = w(t) => F(t,O)
O(to)=O0
mi sono scritto l'algoritmo di integrazione (un poo vedendo qui u po altrove)
ora...COME GLIELA PASSO LA FUNZIONE ??
purtroppo nn sono molto bravo in matlab
questo è il codice di runge-kutta
function [t,y] = RK(odefun,tspan,y0,h)
%RK Metodo di Runge?Kutta IV per ODE
% SINTASSI
% [t,y] = RK(odefun,tspan,y0,h)
% INPUT
% odefun : La funzione f(t,y) del problema differenzale
% y'= f(t,y)
% tspan : Valore iniziale tspan(1) e finale tspan(2)
% della procedura di integrazione
% y0 : Condizione initiale var. dipendente
% h : passo temporale
% OUTPUT
% t : Vettore dei valori delle variabili indipendenti
% y : Vettore dei valori delle variabili dipendenti
% pre?allocazione memoria var. indipendenti
t = (tspan(1):h:tspan(2))';
% numero totale di elementi
n = length(t);
% pre?allocazione memoria var. dipendenti
y = zeros(n,1);
% valore iniziale
y(1) = y0;
% ciclo principale
for j=1:n-1
k1 = feval(odefun,t(j),y(j));
k2 = feval(odefun,t(j)+h/2,y(j)+k1*(h/2));
k3 = feval(odefun,t(j)+h/2,y(j)+k2*(h/2));
k4 = feval(odefun,t(j)+h,y(j)+h*k3);
y(j+1) = y(j)+(h/6)*(k1+2*k2+2*k3+k4);
end
vi prego aiuto!!!
io devo effettuare l'integrazione con runge-kutta
dO/dt = w(t) => F(t,O)
O(to)=O0
mi sono scritto l'algoritmo di integrazione (un poo vedendo qui u po altrove)
ora...COME GLIELA PASSO LA FUNZIONE ??
purtroppo nn sono molto bravo in matlab
questo è il codice di runge-kutta
function [t,y] = RK(odefun,tspan,y0,h)
%RK Metodo di Runge?Kutta IV per ODE
% SINTASSI
% [t,y] = RK(odefun,tspan,y0,h)
% INPUT
% odefun : La funzione f(t,y) del problema differenzale
% y'= f(t,y)
% tspan : Valore iniziale tspan(1) e finale tspan(2)
% della procedura di integrazione
% y0 : Condizione initiale var. dipendente
% h : passo temporale
% OUTPUT
% t : Vettore dei valori delle variabili indipendenti
% y : Vettore dei valori delle variabili dipendenti
% pre?allocazione memoria var. indipendenti
t = (tspan(1):h:tspan(2))';
% numero totale di elementi
n = length(t);
% pre?allocazione memoria var. dipendenti
y = zeros(n,1);
% valore iniziale
y(1) = y0;
% ciclo principale
for j=1:n-1
k1 = feval(odefun,t(j),y(j));
k2 = feval(odefun,t(j)+h/2,y(j)+k1*(h/2));
k3 = feval(odefun,t(j)+h/2,y(j)+k2*(h/2));
k4 = feval(odefun,t(j)+h,y(j)+h*k3);
y(j+1) = y(j)+(h/6)*(k1+2*k2+2*k3+k4);
end
vi prego aiuto!!!
Devi definire una function che rappresenta la funzione della tua equazione differenziale, quindi nel tuo caso $f(t,y(t)) = w(t)$ e quindi in matlab dovrai definire qualcosa del tipo
$f=\text{@}(t,y) w(t)$, dove $w$ è la funzione che hai nel problema.
Poi se devi creare il vettore $tspan$ che contiene come primo elemento il primo istante temporale $t_0$ e l'ultimo, quindi se vuoi la soluzione su $[t_0 t_n]$, $tspan=[t0 tn]$ dove $t0,tn$ sono dei valori scalari che indicano gli istanti temporali suddetti.
A questo punto puoi usare la tua function RK.
$f=\text{@}(t,y) w(t)$, dove $w$ è la funzione che hai nel problema.
Poi se devi creare il vettore $tspan$ che contiene come primo elemento il primo istante temporale $t_0$ e l'ultimo, quindi se vuoi la soluzione su $[t_0 t_n]$, $tspan=[t0 tn]$ dove $t0,tn$ sono dei valori scalari che indicano gli istanti temporali suddetti.
A questo punto puoi usare la tua function RK.
grazie per la risp...ti dico pero il problema....io w(t) non la conosco analiticamente...so so i valori che assume in determinati istanti di tempo...come facciio??:(((
puoi provare ad interpolare i dati che hai e usare l'interpolatore come $w(t)$ nell'equazione differenziale, tipo usare l'interpolazione spline che garantisce che la funzione interpolante sia di classe $C^2$.
e se invece gli passassi i i valori negli istanti di tempo che ho io?
Il metodo RK come vedi valuta la $f(t,y)$ non solo negli istanti $t_k=t_0 + k\cdot h$, ma anche a valori intermedi, dunque dovresti avere tutti questi punti, non è semplice la cosa.
quindi mi consigli di fare l'interpolazione e usare la funzione che ottengo. giusto?
grazie mille x il consiglio
grazie mille x il consiglio
Considera che comunque si sta introducendo un errore aggiunto causato dall'interpolazione dei dati, però non credo si possa fare molto altro, anche perchè se l'obiettivo è trovare una soluzione al problema di cauchy, ecco considera che la soluzione che calcoli ha un qualche senso nell'intervalo in cui sono contenuti i dati, fuori da quell'intervallo si perde completamente di significato.
Se i dati che hai sono abbastanza fitti e quasi equidistanti potresti anche provare a costruire una sorta di eulero esplicito sugli intervalli data dai dati che possiedi, alla fine la tua equazione differenziale si riduce al calcolo di una primitiva di $w(t)$.
Se i dati che hai sono abbastanza fitti e quasi equidistanti potresti anche provare a costruire una sorta di eulero esplicito sugli intervalli data dai dati che possiedi, alla fine la tua equazione differenziale si riduce al calcolo di una primitiva di $w(t)$.
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo