Soluzione approssimata, equazione differenziale runge-kutta

[_johnnyfreak_]

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:
\( \displaystyle \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}} \)

L'algoritmo che sto studiando è RK4:

Per risolvere \( \displaystyle \frac{dy}{dt} = f(t, x) \)
\( \displaystyle h = \frac{x_{finale} - x_{iniziale}}{N} \) è il passo di integrazione
\( \displaystyle k_1 = f(x_n, y_n) \)
\( \displaystyle k_2 = f(x_n + \frac{h}{2}, y_n +\frac{k_1*h}{2}) \)
\( \displaystyle k_3 = f(x_n + \frac{h}{2}, y_n + \frac{k_2*h}{2}) \)
\( \displaystyle k_4 = f(x_n + h, y_n + k_3*h) \)
\( \displaystyle y_{n+1} = y_n + \frac{k_1}{6} + \frac{k_2}{3} + \frac{k_3}{3} + \frac{k_4}{6} + O(h^5) \)

Ma non riesco ad applicarlo perchè \( \displaystyle {f{{2}}} \) ha più variabili.

Potreste darmi una mano?
Grazie

[Ska]

Stai lavorando con un sistema di equazioni differenziali del primo ordine, dunque nella forma

\( \displaystyle \begin{cases} \underline y'(t) = \underline F(t, \underline y(t)) \\ \underline y(t_0) = \underline y_0\end{cases} \)

Nel tuo caso \( \displaystyle \underline y(t) = \begin{bmatrix}\theta(t) \\ \omega(t)\end{bmatrix} \) e \( \displaystyle F(t,\underline y(t)) = \begin{bmatrix}\omega(t) \\ A*\sin(\theta(t)) + B*v_a - B*\omega(t)\end{bmatrix} \)

Dunque rispetto al problema generale avrai

\( \displaystyle \underline k_1=\underline F(x_n, \underline y_n) \)
\( \displaystyle \underline k_2=\underline F(x_n + \frac{h}{2}, \underline y_n + \underline k_1 \frac{h}{2}) \)
\( \displaystyle \underline k_3=\underline F(x_n + \frac{h}{2}, \underline y_n + \underline k_2 \frac{h}{2}) \)
\( \displaystyle \underline k_4=\underline F(x_n + h, \underline y_n + \underline k_3 h) \)
\( \displaystyle 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) \)

[_johnnyfreak_]

Grazie per la risposta.
Se ho capito bene devo iterare il processo per ogni funzione (quindi calcolare 2 \( \displaystyle {k}_{{i}} \) e 2 \( \displaystyle {y}_{{{i}+{1}}} \) ad ogni passo.

Ho provato a scrivere il codice in matlab:

Questa è la funzione che applica l'algoritmo:

Testo nascosto, fai click qui per vederlo

Codice: Seleziona tutto


function [ ] = myRunge( f, g, tmin, tmax, N, x0, y0 )

gr = 9.8;
m = 10;
ma = 0.1;
va = 0.1;
r = 0.1;
d = 0.1;

A = 3 * gr / (2 * r)
B = 3 * ma * d / (2 * m * r * r)


h = (tmax - tmin)/N;
t = tmin:h:tmax;

x = zeros(1, length(t));
y = zeros(1, length(t));

asin_Arg = B*va/A;

if(asin_Arg > 1 || asin_Arg < -1)
  ME = MException('VerifyOutput:OutOfBounds', 'errore :) ');
   throw(ME);
end

x(1) = asin(asin_Arg);%x0;
y(1) = y0;


for index = 1:(length(t)-1)
    xi = x(index);
    yi = y(index);
   
    k11 = h*feval(g, xi, yi);
    k22 = h*feval(g, xi+h/2, yi+k11/2);
    k33 = h*feval(g, xi + h/2, yi + k22/2);
    k44 = h*feval(g, xi + h, yi + k33);
   
    k1 = h*feval(f, xi, yi, A, B, va);
    k2 = h*feval(f, xi+h/2, yi+k1/2, A, B, va);
    k3 = h*feval(f, xi + h/2, yi + k2/2, A, B, va);
    k4 = h*feval(f, xi + h, yi + k3, A, B, va);

    x(index+1) = xi+ (k11 + 2*k22 + 2*k33 + k44)/6;
    y(index+1) = yi + (k1 + 2*k2 + 2*k3 + k4)/6;
end

hold off
plot(t, y);%accel
hold on
plot(t, x, 'r');%velocità angolare
plot(x, y, 'g');
end


Queste due sono le funzioni del sistema:

Testo nascosto, fai click qui per vederlo

Codice: Seleziona tutto

function [ y ] = mffunc( x, w , A, B, va)
y = -A*sin(x) + B*va - B*w;
end

function [z] = g(x, y)
    z = y;
end


Che chiamo in questo modo (il penultimo parametro non serve e devo eliminarlo):

Codice: Seleziona tutto

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?

[Ska]

No no, puoi fare tutto manipolando direttamente i vettori

Codice: Seleziona tutto

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 \( \displaystyle t \) e \( \displaystyle \underline y \)
- 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 \( \displaystyle {u} \) di dimensione \( \displaystyle {n}\times{T} \), dove \( \displaystyle {n} \) è il numero di componenti della variabile \( \displaystyle {y} \) e \( \displaystyle {T} \) è il numero degli istanti temporali considerati.

Si inizializza la prima colonna della matrice con \( \displaystyle {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 \( \displaystyle {u} \) per mantenere la convenzione che il risultato contiene nelle colonne le soluzioni di ciascuna componente.

La tua funzione \( \displaystyle {F} \) da passare dovrebbe essere ad esempio

Codice: Seleziona tutto

F=@(t,y)[y(2); A*sin(y(1))+B*v_a-B*y(2)];


[_johnnyfreak_]

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?

[Ska]

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

[_johnnyfreak_]

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.

[oniudra]

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)[/b][/b]
%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!!!

[oniudra]

la mia mail è ardu.cicchetti@gmail.com

vi prego datemi una mano!!
grazie a tutti!!

[Ska]

Devi definire una function che rappresenta la funzione della tua equazione differenziale, quindi nel tuo caso \( \displaystyle {f{{\left({t},{y}{\left({t}\right)}\right)}}}={w}{\left({t}\right)} \) e quindi in matlab dovrai definire qualcosa del tipo

\( \displaystyle {f{=}}\text{@}{\left({t},{y}\right)}{w}{\left({t}\right)} \), dove \( \displaystyle {w} \) è la funzione che hai nel problema.

Poi se devi creare il vettore \( \displaystyle {t}{s}{p}{a}{n} \) che contiene come primo elemento il primo istante temporale \( \displaystyle {t}_{{0}} \) e l'ultimo, quindi se vuoi la soluzione su \( \displaystyle {\left[{t}_{{0}}{t}_{{n}}\right]} \), \( \displaystyle {t}{s}{p}{a}{n}={\left[{t}{0}{t}{n}\right]} \) dove \( \displaystyle {t}{0},{t}{n} \) sono dei valori scalari che indicano gli istanti temporali suddetti.

A questo punto puoi usare la tua function RK.