Laboratorio 9

Equazioni differenziali ordinarie

Per risolvere numericamente un problema di Cauchy per un sistema di equazioni differenziali ordinarie del primo ordine con condizioni iniziali assegnate - dopo aver controllato che ci siano le condizioni teoriche di esistenza ed unicità in piccolo ! - esistono innumerevoli algoritmi, ciascuno adatto a particolari classi di problemi, e sicuramente molti altri algoritmi aspettano ancora di essere inventati.

Ricordiamo comunque i metodi elementari come quello di Eulero, e quelli grafici come il metodo delle isocline, che si implementano senza sforzo.(vedi questi richiami)

La scelta dell'algoritmo va fatta oculatamente e non ci sono ricette magiche universali.
Alcuni algoritmi usano un passo fisso, ma la maggior parte usa un passo variabile adattivo.

Molti di questi algoritmi sono disponibili sia in octave che matlab ma le implementazioni e il modo di usarli possono essere diversi.

In octave una routine general purpose è comunque lsode, che con il comando di configurazione lsode_options può essere adattata a diverse situazioni comuni, anche variando l'algoritmo di quadratura.

Per esempio consideriamo il problema di Cauchy:

x'(t)=(3-x(t))* x(t) + sin(t)
x(0)=1

Per calcolare la soluzione per x in [0,5] procediamo in questo modo:
f=@(x,t) (3-x).*x + sin(t) tr=linspace(0,5,500); xr=lsode(f,1,tr); tl=linspace(0,-5,500); xl=lsode(f,1,tl); plot(tl,xl,tr,xr);

Si noti che i valori di t devono iniziare con il punto iniziale in cui viene assegnato il valore di x.
Se si vuole calcolare la soluzione in un intervallo che contiene al suo interno l'istante iniziale si deve fare due distinte operazioni, una a destra e l'altra a sinistra del punto iniziale.

In matlab (e anche in octave) per lo stesso problema di Cauchy si può usare per esempio ode45:
f=@(t,x) (3-x).*x + sin(t) [tr,xr]=ode45(f,[0,5],1); [tl,xl]=ode45(f,[0,-5],1); dt=tr[2:end]-tr[1:end-1]; plot(tl,xl,tr,xr);

Si noti che:

L'ordine in cui sono messe le variabile nelle funzione deve essere invertito !!!!
ode45 crea i suoi passi che spesso sono di lunghezze variabili (vedi dt ...)
se non si assegna il valore di output di ode45, viene automaticamente fatto il plot
per cambiare i parametri si può usare odeset
nelle versioni più recenti di octave per usare ode45 occorre esplicitamente caricare il pacchetto con il comando pkg load odepkg

Per esempio:
f=@(t,x) (3-x).*x + sin(t) options=odeset('MaxStep',0.01,'InitialStep',0.01); ode45(f,[0,20],1,options);

Infine in scilab una routine general purpose è ode e anche qui la sintassi è lievemente diversa ma simile:
deff('u=fz(t,x)','u=(3-x).*x + sin(t)'); tr=linspace(0,20,500); xr=ode(1,0,tr,fz); plot(tr,xr);

Si noti il modo diverso di definire le funzioni inline, e la sintassi di ode che ha come argomenti il valore iniziale di x, di t, e poi come in lsode il vettore delle variabili indipendenti e la funzione - che però prende gli argomenti nell'ordine di ode45.

Laboratorio Didattico di Matematica Computazionale - Sergio Steffè - AA 2017/2018 - PISA