Metodi numerici per equazioni differenziali ordinarie - laboratorio
Docente: Beatrice Meini

Lezione 2: I metodi di Runge-Kutta

• La function function [t, x] = RK(fname, tspan, x0, h, G, beta) applica un metodo di Runge-Kutta esplicito per risolvere un IVP. La variabile fname rappresenta la funzione f(t,x) che definisce l'equazione differenziale, tspan un vettore di dimensione 2 che rappresenta l'intervallo d'integrazione, x0 e' il valore iniziale, h e' il passo di discretizzazione, G e beta sono la matrice Gamma e il vettore Beta che definiscono il metodo di RK. In uscita, t contiene l'intervallo discretizzato, e x le approssimazioni della soluzione.
  La function RK.m utilizza la function [tnew, xnew, fnew] = RKstep(fname, t, x, fx, h, G, ro, beta) che applica un passo del metodo di Runge-Kutta, definito dai parametri G, ro e beta.
  Il file G_beta.tar contiene le matrici Gamma e i vettori beta che definiscono vari metodi di RK: il metodo di Eulero, il metodo di Heun, un metodo a 3 stadi, e il classico metodo di RK a 4 stadi. Salvare il file e scomprimere l'archivio con il comando tar xvf G_beta.tar. Otterrete dei file .mat (ad esempio GHeun.mat, bHeun.mat), che caricate in Matlab con il comando load [nome].mat. I metodi definiti da questi parametri sono tali che un metodo a m stadi ha ordine di converegnza m.

• La function [y] = f3(t,x) e' la funzione che definisce l'equazione differenziale x'(t)= -x+2 x^2 exp(-t). La soluzione dell'IVP con condizione iniziale x(0)=1/3, nell'intervallo [0 T], e' la funzione definita dalla function [y] = sol3(t). Approssimare la soluzione con i vari metodi di RK, disegnare la soluzione e calcolare l'errore, ad esempio con i comandi:
  h=1.e-2;
  [t1,x1]=RK(@f3, [0 1], 1/3, h, GEulero, bEulero);
  [t2,x2]=RK(@f3, [0 1], 1/3, h, GHeun, bHeun);
  [t3,x3]=RK(@f3, [0 1], 1/3, h, GRK3, bRK3);
  [t4,x4]=RK(@f3, [0 1], 1/3, h, GRK4, bRK4);
  y=sol3(t1);
  plot(t1,y,t1,x1,t2,x2,t3,x3,t4,x4)
  legend('esatta','Eul','Heun','RK3','RK4')
  err1=abs(x1-y);
  err2=abs(x2-y);
  err3=abs(x3-y);
  err4=abs(x4-y);
  semilogy(t1,err1,t2,err2,t3,err3,t4,err4)
  legend('Eul','Heun','RK3','RK4')

• Costruire una function che definisce l'equazione differenziale x'(t)=-200 (x(t)+ sin(t)). Risolvere l'IVP con condizione iniziale x(0)=2, nell'intervallo [0 5], mediante i 4 metodi di Runge-Kutta, con passo h=1.e-2. Disegnare, sulla stesso grafico, le approssimazioni ottenute. Che cosa si osserva? Perche'? Modificare il passo di discretizzazione e fare altre prove.

• La function [y] = moto(t,x) e' la funzione che definisce una equazione differenziale che descrive la traiettoria di un corpo di massa trascurabile, soggetto all'attrazione di altri 2 corpi, uno di massa mu (piccola) e uno di massa 1-mu, che si muovono in un'orbita planare. L'equazione differenziale e' descritta nell'Esercizio 4.12 del libro di Ascher, Petzold, che trovate qui. Approssimare la soluzione con il metodo di RK a 4 stadi, prendendo come condizione iniziale x0=[0.994; 0; 0; -2.00158510], e intervallo [0, 17.1]. Tracciare il grafico della traiettoria del corpo:
  
  [t,x]=RK(@moto, [0, 17.1], x0, 1.e-2, GRK4, bRK4);
  plot(x(:,1),x(:,2))
  
  Cosa si osserva? Provare a diminuire il passo di discretizzazione.
  Risolvere lo stesso problema con il risolutore ode45 di Matlab (si dia il comando help ode45):
  
  [t1,x1]=ode45(@moto, [0, 17.1], x0)
  plot(x(:,1),x(:,2),'+')
  Quanti, e dove sono maggiormente concentrati, i punti di discretizzazione dell'intervallo?