METODO KUTTA (2°ordine)

Spazio studenti
1

Lo so che è un’impresa disperata ma se qualcuno sa aiutarmi gli faccio una statua!

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)=Acos(wt+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’’+(w^2)y=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_analisi/pdf/rungekutta.pdf

2

Uhmm ora non ho Matlab dove sono per provare… :kissing_heart:
L’unico consiglio che ti posso dare è di rivedere intervalli, campionamenti, frequenze e metodi di calcolo… se non lo hai già fatto… prova a cambiare e provare un po’ per volta e vedere cosa succede… mi è capitato spesso che simulazioni andassero a farsi benedire per errori numerici e non procedurali, però la risoluzione era quasi sempre “unica” caso per caso… :?
Altrimenti, aspetta il contributo di altri utenti… credo che il problema non sia così “anomalo”…

3

grazi della risposta ma è tutta la mattina che prova a cambiare step di calcolo frequenze e anche le condizioni iniziali ma il risultato è sempre incoerente!

La conclusione a cui sono arrivato è che o non ho ben capito questo metodo numerico sbagliando nell’impostazione, oppure la precisione di questo metodo è del millesimo di secondo (cosa che un po mi sconcerta)!!

Non è cosi impreciso questo metodo vero?Anche perche se uso le funzioni di matlab tipo ODE45 che non è altro che questo metodo ma svilupato fino al 4° ordine, i risultati analitici/numerici sono pochissimo differenti qualsiasi intervallo temporale prendo (anche 50secondi)!!! :incazzato: :incazzato:

4

Di quanto diverge? Diversi ordini di grandezza?

Suggerimento: modifica il valore iniziale e il passo da 0.01 a 0.02.

Motivazione: 0.01 non ha una rappresentazione finita in una codifica binaria di numeri decimali a virgola mobile, che presumo venga impiegata anche da Matlab. Le imprecisioni nella rappresentazione possono dunque propagarsi più rapidamente.

Paolo Amoroso

5

in 10 secondi mi sembra che diverge di 10 alla 300 :incazzato: :incazzato:

cmq ho provato a risolverla con ode23 la funzione intrinseca di matlab che nella descrizione dice essere Runge-Kutta al secondo ordine (ossia quello che sto cercando di fare)…indovinate???La soluzione è perfetta,sbaglia di un millessimo, quindi c’è qualcosa che non va alla base della mia impostazione!

6

Credo di aver trovato il problema: in pratica tu in queste righe
…
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);
…
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;
…
moltiplichi per t, ma t non è il passo di integrazione, ma il tempo proprio (che quindi incrementa sempre di volta in volta). Credo (da quel che ora ricordo sui metodi numerici) che debba usare h=0.01 (nel tuo caso) invece di t.
Prova e fammi sapere!!!

7

Purtroppo non ho tempo in questi giorni, ma nel gruppo delle ODE implementate in matlab usa la ode113 che è una delle migliori

8

AJ lo scopo è proprio quello di non usare i solver di matlab (sarebbe troppo semplice)!!

ArTax scusa ma non ho capito bene, secondo te dov’è l’errore?Cioè come ho fatto io invece di step di tempo uguali mano mano viene aumentato??

9

Ah, ok, sorry

10

% 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;

% 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:10
i=i+1;
h=0.01;
% vettore con gli istanti di tempo
T(i+1,1)=t;
% passo intermedio dopo ho imposto t/2
mx(i+1,1)=mx(i,1)+my(i,1)*(h/2);
my(i+1,1)=my(i,1)-(w^2)mx(i,1)(h/2);

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

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

end

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

Ho capito cosa intendevi con il passo e ho corretto in questo modo!Adesso viene ma non è proprio identico, mi viene comunque una sinusoide ma dopo un po di tempo (ma adesso posso integrare anche per 10 secondi) i massimi delle onde si sfasano proprio!

Dici che questo è una conseguenza dell’imprecisione del calcolo o c’è ancora qualche errore di programmazione?

11

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

Questo è il plottaggio

12

% 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)=Acos(wt+phi);

end

plot(T,x,T,z);

ECCOLOOOOOO Adesso mi viene posto un grafico!

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

GRAZIE a tutti, ora praticamente devo iniziare a fare il mio elaborato per il sistema vibrante! :shocked:

13

Quando avrai un po’ più di tempo ti suggerisco la lettura del capitolo sul debugging del libro The Practice of Programming di Brian Kernighan e Rob Pike. Quei consigli sono la sintesi della saggezza dei programmatori Jedi.

Paolo Amoroso