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

[valentino861]

Si mi serve per meccanica applicata alle macchine la parte sulle vibrazioni!!


% DATI:

m=2;
k=2;
c=0.5;
omega=0.8;
F=1;

% CONDIZIONI INIZIALI:

xo=0;
vo=10;
x(1,1)=xo;
y(1,1)=vo;
T(1,1)=0;
i=0;
% COSTANTI:

% pulsazione naturale
omegan=sqrt(k/m);
% coefficiente smorzamento
S=c/(2*m*omegan);
% pulasazione smorzata
omegad=omegan*sqrt(1-S^2);
% ampiezza
X=sqrt(xo^2+((vo+S*xo*omegan)/omegad)^2)
% fase iniziale
phi=atan((vo+S*xo*omegan)/(xo*omegad))
% fase iniziale della soluzione particolare
PHI=atan((omega*c)/(k-m*omega^2));
% passo d'integrazione
h=0.01;
for t=0.1:0.01:80
i=i+1;
T(i+1,1)=t;
% SOLUZIONE GENERALE
G(i+1,1)=X*(expm(-S*omegan*t))*cos(omegan*t-phi);
% SOLUZIONE PARTICOLARE
P(i+1,1)=(F*sin(omega*t-PHI))/(sqrt((k-m*omega^2)^2+(omega*c)^2));
% SOLUZIONE NUMERICA CON IL METODO DI KUTTA (2°ordine)

% valori mx e my sono quelli calcolati con h/2:
mx(i+1,1)=x(i,1)+y(i,1)*(h/2);
my(i+1,1)=y(i,1)+((F/m*sin(omega*t))-(2*S*omegan*y(i,1))-(omegan^2*x(i,1)))*(h/2);
% valori x e y sono quelli per h:
x(i+1,1)=x(i,1)+my(i+1,1)*h;
y(i+1,1)=y(i,1)+((F/m*sin(omega*(t+h/2)))-(2*S*omegan*my(i+1,1))-(omegan^2*mx(i+1,1)))*h;
end

% SOLUZIONE ANALITICA ESATTA
sol=P+G;

% GRAFICO
plot(T,x,T,sol);


Lo so in alcuni punti è leggermente ridondante ma per la parte analitica non posso farci niente, mi sono risolto da solo l'equazione differenziale per quelle costanti e condizioni iniziali e poi l'ho semplicemente scritta.

Mentre per il metodo di Kutta oggi il prof mi ha detto che devo fare per forza cosi, ossia calcolare l'incremento con uno step intermedio (h/2),e poi calcolare il nuovo valore reale con l'incremento cosi ottenuto!

Posto il grafico di come lo consegnerò domani al prof per presentarmi all'orale

ps:oggi ho fatto lo scritto...speriamo bene!!!!

[valentino861]

Equazione differenziale dell'oscillatore

$mx''+cx'+kx=Fcos(\omegat)$

Ho diviso l'equazione del secondo ordine in due del primo ordine

$x'=y$
$y'=-(c/m)y-(k/m)x+F/mcos(\omegat)$

Primo sistema per calcolare l'incremento

$x(t+h/2)= x(t)+y(t)h/2$
$y(t+h/2)= y(t)+(-(c/m)y(t)-(k/m)x(t)+F/mcos(\omegat))h/2$

Sistema per calcolare la funzione:

$x(t+h)=x(t)+y(t+h/2)h$
$y(t+h)=y(t)+(-(c/m)y(t+h/2)-(k/m)x(t+h/2)+F/mcos(\omega(t+h/2)))h$