METODO KUTTA (2° ORDINE)

[valentino861]

Sono disperato....devo consegnare un elaborato e non ci capisco niente con questo metodo numerico chi mi aiuta??

Devo integrare un'equazione differenziale per un sistema vibrante (ma adesso per capire come funzione ci possiamo anche basare sull'equazione di un pendolo semplice senza smorzamento)!

Il problema è che la soluzione numerica diverge tantissimo dalla soluzione analitica e non so quanto sia giusto!Cioè se prendo degli step di calcolo di un centesimo di secondo...gia dopo mezzo secondo i risultati sono allucinanti.

Vi posto il programma che sto cercando di fare con matlab ma il problema credo si a monte!


Allora io ho usato il metodo di Eulero per trovarmi il valore della funzione dopo la metà del tempo,per capirci:

% sono le condizioni iniziali
mx(1,1)=0.5;
my(1,1)=0;
x(1,1)=mx(1,1);
y(1,1)=my(1,1);

% pulsazione del pendolo
w=4;
i=0;
% costanti che mi sono trovato per la soluzione analitica
A=x(1,1);
phi=0;
z(1,1)=0.5;


for t=0.01:0.01:0.2
i=i+1;
% vettore con gli istanti di tempo
T(i+1,1)=t;
% passo intermedio ho imposto t/2
mx(i+1,1)=mx(i,1)+my(i,1)*(t/2);
my(i+1,1)=my(i,1)-(w^2)*mx(i,1)*(t/2);

% questo dovrebbe essere l'integrazione secondo Kutta
x(i+1,1)=x(i,1)+my(i+1,1)*t;
y(i+1,1)=y(i,1)-(w^2)*mx(i+1,1)*t;
% soluzione analitica

z(i+1,1)=A*cos(w*t+phi);

end

T(1,1)=0;
% soluzione analitica
subplot(2,1,1);
plot(T,z);
% soluzione numerica (Kutta)
subplot(2,1,2);
plot(T,x);



Per chi ha matlab basta prenderlo e copiarlo anche sul command window.

Ovviamente sto parlando di:

$y''+\omega^2y=0$




Allora che ve ne pare dove ho sbagliato?La soluzione analitica è giusta perche se faccio un plottaggio in un arco di tempo più ampio mi viene una sinusoide, mentre se plotto la soluzione numerica anche per un tempo di un secondo mi diverge A BESTIA!!



PS:dopo il metodo di Eulero in questo link è accennato a ciò di cui sto parlando

http://www.afs.enea.it/gianness/corso_a ... ekutta.pdf

[amel3]

Hai fatto un po' di pasticci a dir la verità, non per scoraggiarti, ma mi sembra un tantino ridondante tra l'altro...
Comunque l'errore principale credo (credo, eh) che sia nel fatto che hai scelto un passo variabile invece che un passo fisso (che "classicamente" si indica con $h$).
Se devi fare una cosa al volo, non posso aiutarti, mi spiace, ma se hai qualche giorno (così impari di più e con più calma) nei limiti delle mie possibilità posso darti una mano.

[valentino861]

ciao e grazie della risposta!Potresti darmi qualche dritta allora?non capisco cosa dici sullo step di calcolo!

[amel3]

% sono le condizioni iniziali
mx(1,1)=0.5;
my(1,1)=0;
x(1,1)=mx(1,1);
y(1,1)=my(1,1);

% pulsazione del pendolo
w=4;
i=0;
% costanti che mi sono trovato per la soluzione analitica
A=x(1,1);
phi=0;
z(1,1)=0.5;

%passo integrazione
h=0.01;

for t=0.01:h:0.2
i=i+1;

% questo dovrebbe essere l'integrazione secondo Kutta
x(i+1,1)=x(i,1)+my(i,1)*h;
y(i+1,1)=y(i,1)-(w^2)*mx(i,1)*h;
% soluzione analitica

z(i+1,1)=A*cos(w*t(i+1)+phi);

end

% soluzione analitica
subplot(2,1,1);
plot(t,z);
% soluzione numerica (Kutta)
subplot(2,1,2);
plot(t,x); 

Prova a vedere intanto se così funzia (non è detto, l'ho corretto al volo...), in questo momento non posso usare Matlab.

Però io lo riscriverei quasi tutto 'sto script sinceramente....

[valentino861]

Infatti l'ho corretto in modo simile e ora mi viene però a lungo andare mi continua a divergere è normale?

http://img370.imageshack.us/img370/8997/45209323zf2.jpg

Questo è un plottaggio per $\omega=1$ e in 20 secondi di tempo....è normale quasto andamento?

[amel3]

confrontali insieme, così:

plot(t,z,t,x);

[amel3]

A pensarci bene ho dimenticato un pezzo, come ha fatto a venirti?

Ora ho corretto...

[valentino861]

avevo gia fatto la correzione un secondo prima che tu postassi la tua!!

In ogni caso nella tua versione c'è un t(i,1) questo vettore non è mai stato inizializzato i ci ho messo direttamente t no??

[Conte_De_Saint_venant]

%Funzione che risolve un problema ai valori iniziali con il metodo di
%rungekutta esplicito
%[y,t]=eulero(a,b,f,yo)
%a,b =estremi dell'intervallo
%f funzione da valutare f(t,y'), inserire la funzione tra apici
%y soluzione calcolata sul vettore t delle ascisse
%intervallo su cui viene calcolata la soluzione

function [y,t]=rungekutta(a,b,f,yo)
%
%calcolo dei parametri necessari per l'analisi
clc;
discr=input('Inserire in quante parti discretizzare l''intervallo: ');
t=linspace(a,b,discr);
h=t(2)-t(1);
y(1)=yo;
dt=t(2)-t(1);
%
fun=inline(f,'t','y');

%funzione che calcola la soluzione
for j=1:1:length(t)-1
K1=fun(t(j),y(j));
K2=fun((t(j)+h/2),(y(j)+h/2*K1));
K3=fun((t(j)+h/2),(y(j)+h/2*K2));
K4=fun((t(j)+h),(y(j)+h*K3))
y(j+1)=y(j)+dt/6*(K1+2*K2+2*K3+K4);
end;






Guarda questo listato l'avevo scritto tempo fa per il metodo del quarto ordine. se hai sotto mano i tableau di butcher lo puoi modificare facilmente
per il metodo del secondo ordine.

Studi ingegneria? Per la risposta di sistemi vibranti comunque caricati è efficiente il metodo di Newmark(1959) ....

[valentino861]

questa è la soluzione analitica e la soluzione che hai postato tu del ciclo for!

Ma la tua versione non è la normale soluzione di Eulero??Non presenta il passo intermedio caratteristiche del metodo di Kutta!

[amel3]

Cioè hai fatto z=A...eccetera ecc. tirandolo fuori dal for? In tal caso ok.

[valentino861]

% sono le condizioni iniziali
x(1,1)=0.5;
y(1,1)=0;
% pulsazione del pendolo
w=5;
% costanti che mi sono trovato per la soluzione analitica
A=x(1,1);
phi=0;
z(1,1)=0.5;

i=0;
T(1,1)=0;
h=0.01;
for t=0.01:0.01:5
i=i+1;
% vettore con gli istanti di tempo
T(i+1,1)=t;
% questo dovrebbe essere l'integrazione secondo Kutta
x(i+1,1)=x(i,1)+[y(i,1)-(w^2)*x(i,1)*(h/2)]*h;
y(i+1,1)=y(i,1)-(w^2)*x(i+1,1)*h;
% soluzione analitica

z(i+1,1)=A*cos(w*t+phi);

end


plot(T,x,T,z);



ECCOLOOOOOO Adesso mi viene posto un grafico!

http://img84.imageshack.us/img84/382/mat2ir0.jpg

[amel3]

Strano però i grafici vengano così diversi; prova a raffinare il passo. Vabbè che il metodo converge alla soluzione esatta solo con ordine 1, però...

(P.S.: Se provi Runge-Kutta del IV ordine dovresti ottenere una soluzione quasi indistinguibile da quella esatta; questo non c'entra con quello di cui stiamo parlando ora ovviamente)

[Conte_De_Saint_venant]

Non ho capito tanto che mi hai scritto....cmq il mio è un metodo di Runge e Kutta del 4° ordine....Il metodo di eulero implicito ( meglio conosciuto come metodo di Heun) ed il metodo di eulero esplicito sono due casi particolari degli infiniti metodi di Runge e Kutta a due stadi espliciti.....
Detto questo non mi sembra che il tuo algoritmo possa funzionare non do un senso hai ruoli di my....
Cmq guarda:

u(i+1)=u(i)+ h (b1*K1 + b2 *K2)
K1=f(t(i),u(i));
K2=f(t(i)+ h c2, y(i) + h c2 K1);

deve essere b1+b2=1, c2*b2=1/2..... scelgo ad esempio b1=1/2; b2=1/2; c2=1;

allora il ciclo è questo:
%fun è la funzione da valutare y'= fun(t,y);
% t è un vettore equispaziato di h (es t=linspace (0,30,300); h=t(2)-t(1))
%u(1)= la tua condizione iniziale...
for i=1:1:length(t)-1

K1= fun(t(i),u(i));

K2=fun(t(i)+ h , u(i) + h* K1);

u(i+1)=u(i)+ h/2* (K1 + K2);

end

[amel3]

Scusa mi sa che abbiamo fatto un po' di confusione... io parlavo con valentino86 e dormivo abbastanza da non aver notato per nulla che avevi postato anche tu (sic!!). Comunque occhio che Heun ha ordine 2 e non 4.

P.S.: Che verso fa un maiale quando tenta di risolvere una edo? Heunk, Heunk!
OK, basta...

[Conte_De_Saint_venant]

Io mi riferivo a Valentino86....ieri ho postato il metodo di runge kutta a 4 stadi e ordine di convergenza 4 perchè avevo un matlab's function che avevo scritto tempo fa.... Lui mi aveva risposto in un modo che non ho capito dicendo che quello era il metodo di Eulero.... Per cui la risposta era per lui.
Successivamente quando mi riferivo ai metodi a due stadi non ho parlato affatto di convergenza (che è due) bensi solamente di un metodo per ricavarlo e ho scritto il relativo ciclo for...

[amel3]

Ok, però attenzione che di solito quando si parla di ordine si intende l'ordine di convergenza, quello che dici tu di solito è chiamato numero di passi o di stadi.
Vabbè, scusa per aver fatto il saputello saccente, che tra l'altro non mi si addice...
Ciao passo e chiudo.

[Conte_De_Saint_venant]

boh .....non ho affatto parlato di convergenza e credo di sapere la differenza tra ordine di convergenza e passi (o stadi). Spero di spiegarmi così:

1) Ho postato un metodo di Runge e Kutta multistep a 4 stadi, che ha convergenza 4
2) Ho semplicemente risposto al quesito di valentno86 che chiedeva un metodo a due stadi.

Cmq l'ordine di convergenza a seconda dello stadio è:
s=1 ===> 1
s=2 ===>2
s=3 ===>3
s=4 ===>4
s=5====>4
s=6 ====>5
s=7 ====>6
s=8 ====> 6

Fonte "Matematica Numerica" di Quarteroni Alfio (Forza ALinghi), Sacco, Saleri dal quale ho studiato a suo tempo metodi numeroci.


Valentino 86 studi ingegneria strutturale?

[valentino861]

Arieccome !!

Allora il programma finale che avevo postato io infatti è sbagliato peche alla fine non è altro che uno sviluppo di taylor al 2°ordine invece che il metodo il Kutta!Infatti oggi il prof mi si mangiato!
In ogni caso la soluzione di Conte_de_Saint_Venant è giusta!In serata mi ci metto è posto la soluzione definitiva del problema che sto analizzando...un oscillatore smorzato sollecitato da una forza armonica!


PS: studio ingegneria meccanica, perchè??

[amel3]

"valentino86":oggi il prof mi si mangiato!

Mi dispiace, non ti volevo far sbagliare, d'altronde ti avevo detto che sarei stato in grado di aiutarti in un po' di tempo, altrimenti era meglio lasciar perdere (avevo pur detto che lo script andava rifatto completamente).
Io ti avevo indicato Eulero onestep perchè credevo dovessi solo cominciare a fare delle prove. Avevo evidentemente frainteso. Poi ero troppo stanco e avevo spento il cervello e non avevo neanche guardato bene il tuo script...
Sicuramente il Conte ti saprà aiutare meglio di me. Cerca anche di capire quello che ti scrive, però, senza copiare bovinamente (non ti offendere, eh )

Scusa ancora, ciao (ora mi zittisco veramente...)

[Conte_De_Saint_venant]

Sono contento sia giusta....comunque veramente consiglio anche a te quel testo se devi trattare dei metodi numerici.
Ti ho chiesto cosa studiavi, per avere magari uno scambio di conoscenze, io studio ingegneria strutturale.
Per cosa ti serve il programmino? Meccanica delle vibrazioni?