Problema metodo shooting.
Salve a tutti!
Sto implementando un algoritmo che permetta di calcolare il periodo di un pendolo date come condizioni al contorno l'angolo iniziale ($\theta_{1}$) al tempo $t_{1}$ e l'angolo finale ($\theta_{N}$) al tempo $t_{2}$. L'equazione del pendolo è $\frac{d^2 \theta}{d t^2} = -\frac{g}{L} sin\theta$ e voglio usare il metodo dello shooting con Runge Kutta al quarto ordine per trovare i valori intermedi e procedere con il calcolo del periodo. Premesso che Runge Kutta l'avevo già scritto per un altro programma ed è funzionante, io ho pensato di fare così:
1. scelgo un valore a caso per $\theta'_{1}(t_{1})$ e calcolo RK.
2. scelgo un altro valore a caso per $\theta'_{2}(t_{1})$ e calcolo RK.
3. a questo punto calcolo la nuova $\theta'(t_{1})$ che andrò ad utilizzare come $\theta'(t_{1}) = \theta'_{1}(t_{1}) + \frac{(\theta'_{2}(t_{1})-\theta'_{1}(t_{1}))}{(\theta_{2}(t_{1})-\theta_{1}(t_{1}))}(\theta_{N}(t_{2}) - \theta_{1}(t_{1}))$
4. Entro in un ciclo e utilizzo la $\theta'_{2}(t_{1})$ e $\theta'(t_{1})$ per calcolare RK, successivamente salvo in $theta'_{2}(t_{1})$ il valore di $\theta'(t_{1})$ e quest'ultima prende il valore dato dall'interpolazione che ho riportato in precedenza.
A questo punto ho scelto di far girare il ciclo finchè non converge al valore finale $\theta_{2}(t_{2})$ oppure non raggiunge un limite di iterazioni.
Il problema è che non capisco perchè a seconda dei valori che do inizialmente converge o meno e in particolare quando converge non mi da sempre gli stessi risultati. Non ho ben compreso il funzionamento di questo metodo oppure ho fatto qualche errore nella scrittura del programma? Grazie a tutti!!
Sto implementando un algoritmo che permetta di calcolare il periodo di un pendolo date come condizioni al contorno l'angolo iniziale ($\theta_{1}$) al tempo $t_{1}$ e l'angolo finale ($\theta_{N}$) al tempo $t_{2}$. L'equazione del pendolo è $\frac{d^2 \theta}{d t^2} = -\frac{g}{L} sin\theta$ e voglio usare il metodo dello shooting con Runge Kutta al quarto ordine per trovare i valori intermedi e procedere con il calcolo del periodo. Premesso che Runge Kutta l'avevo già scritto per un altro programma ed è funzionante, io ho pensato di fare così:
1. scelgo un valore a caso per $\theta'_{1}(t_{1})$ e calcolo RK.
2. scelgo un altro valore a caso per $\theta'_{2}(t_{1})$ e calcolo RK.
3. a questo punto calcolo la nuova $\theta'(t_{1})$ che andrò ad utilizzare come $\theta'(t_{1}) = \theta'_{1}(t_{1}) + \frac{(\theta'_{2}(t_{1})-\theta'_{1}(t_{1}))}{(\theta_{2}(t_{1})-\theta_{1}(t_{1}))}(\theta_{N}(t_{2}) - \theta_{1}(t_{1}))$
4. Entro in un ciclo e utilizzo la $\theta'_{2}(t_{1})$ e $\theta'(t_{1})$ per calcolare RK, successivamente salvo in $theta'_{2}(t_{1})$ il valore di $\theta'(t_{1})$ e quest'ultima prende il valore dato dall'interpolazione che ho riportato in precedenza.
A questo punto ho scelto di far girare il ciclo finchè non converge al valore finale $\theta_{2}(t_{2})$ oppure non raggiunge un limite di iterazioni.
Il problema è che non capisco perchè a seconda dei valori che do inizialmente converge o meno e in particolare quando converge non mi da sempre gli stessi risultati. Non ho ben compreso il funzionamento di questo metodo oppure ho fatto qualche errore nella scrittura del programma? Grazie a tutti!!
Risposte
Guarda, anche io ho lo stesso problema, cerchiamo di risolverlo insieme entro il 16? Altrimenti non lo diamo calcolo..
Con affetto,
Iustin
Con affetto,
Iustin
Ciao, provo ad aiutarti.
Supponiamo di essere sotto l'ipotesi di piccoli angoli, in modo tale da approssimare $sin(\theta) = \theta$.
Supponiamo che l'angolo iniziale sia pari a $\theta_0 = 1° $ e che l'angolo finale sia pari a $\theta_f = 0°$.
Supponiamo ancora di avere un pendolo lungo $0.51 m$; allora il suo periodo sarà, sotto l'ipotesi di piccoli spostamenti, $ T = 1.4s$.
Osserviamo poi che siamo in presenza di un problema ai limiti. Abbiamo, infatti, il valore dell'angolo all'istante iniziale e finale. Ci occorre, per risolvere il problema, il valore della derivata prima all'istante iniziale, per cui dobbiamo implementare una procedura di shooting.
Siccome MatLab gestisce solo equazioni di primo grado, dobbiamo trasformare l'equazione del pendolo in un sistema di due equazioni di primo grado.
Ecco lo script
La funzione riceve in ingresso una matrice theta di due colonne ed una riga, per ogni istante di tempo. Nella prima colonna sono situati tutti i valori della funzione, mentre nella seconda colonna sono conservati tutti i valori della derivata prima della funzione.
Ad esempio, nell'istante di tempo iniziale, la prima colonna conterrà il valore dell'angolo iniziale, mentre la seconda colonna quello della sua derivata prima.
Per gestirlo in MatLab usiamo tuttavia un vettore colonna, ma la sostanza non cambia.
Questa prima function calcola le pendenze, sotto forma di coefficienti K, che servono per il metodo RK di 4 passi.
Creiamo allora lo script che, ricevuto in ingresso il vettore tempo e le condizioni iniziali, implementi il metodo di RK4 utilizzando le K trovate nello script precedente.
Per il metodo di RK4 si può utilizzare il classico array di Butcher
$ \theta_(n+1) = \theta_n + (h)/6*(K1 + 2K2 + 2K3 + K4) $
dove
$ a_(ij) = [ ( 0 , 0, 0, 0),( 0, 1/2, 0, 0),( 0, 0, 1/2, 0),( 0, 0, 0, 1) ] $
$ b = (1/6, 1/3,1/3,1/6) $
$ c = (0, 1/2, 1/2, 1) $
Il metodo di RK4 è implementato da un ciclo for.
La soluzione la mettiamo in una matrice THETA di iTf righe e 2 colonne; come detto sopra, nella prima colonna abbiamo la funzione, mentre nella seconda abbiamo la sua derivata prima (in questo caso spazio e velocità).
Questo script funziona in questo modo;
Il primo istante di tempo è pari alla condizione iniziale; quindi riempio immediatamente la prima riga della matrice THETA con le I.C.
Per gli istanti successivi calcoliamo la funzione THETA, in ogni istante temporale, grazie al metodo RK4.
Supponendo di trovarci all'istante it = 2;
-Prima invio alla function "pendolo" il vettore theta che contiene la condizione iniziale, in modo tale da calcolare la pendenza nel punto iniziale (K1).
-La pendenza nel punto iniziale viene utilizzata per giungere fino al punto $ t + (dt)/2 $;
In questo punto si calcola poi una nuova pendenza (K2), grazie alla function"pendolo", e la si "riporta" di nuovo al punto iniziale
Su questa nuova pendenza K2 ci si sposta fino a raggiungere il punto $ t + (dt)/2 $ dove si calcolerà una nuova pendenza (K3), sempre grazie alla function "pendolo", che si "riporterà" al punto iniziale.
Su quest'altra pendenza K3 ci si sposta fino a raggiungere il punto $ t + dt $ dove si calcolerà una nuova pendenza (K4), ancora grazie alla function "pendolo", pendenza che, insieme alle altre tre precedentemente trovate, verrà utilizzata per aggiornare all'istante temporale successivo il valore della THETA.
Veniamo adesso allo shooting.
Per capire questo script dobbiamo fare un piccolo riepilogo.
Noi stiamo risolvendo l'equazione del secondo ordine, tradotta in un sistema di EDO, con un metodo RK4.
Se avessimo tutte e due le condizioni iniziali, problema risolto!
Invece a noi manca la seconda condizione iniziale che va ricercata con il metodo di shooting.
Assegniamo quindi una certa condizione iniziale, $ s $, sperando che sia giusta. Ma come facciamo a sapere che la condizione iniziale "sparata" sia esatta?
Semplicemente, confrontando il valore della funzione in $ t = iTf $, che viene fuori grazie alla risoluzione con
RK4 più seconda condizione iniziale ricavata con shooting ( THETA(end,1) ), e il valore reale che la funzione
assume in $ t = iTf $ .
Adesso, poiché abbiano salvato le soluzioni del sistema di EDO in una matrice dove ogni i-esima riga corrisponde ad un istante temporale i-esimo ed ogni colonna j-esima rappresenta la rispettiva derivata di ordine j-1, allora diciamo "prendi il valore, contenuto nella matrice delle soluzioni THETA, alla riga ultima e pima colonna e confrontalo col valore che so che la funzione assume all'istante finale".
In questo caso, il valore all'istante finale è nullo ($\theta_f = 0° $), quindi dobbiamo ricercare quel parametro di shooting che rende nullo il valore della funzione THETA in $ t = iTf $.
Per lo shooting useremo un metodo di bisezione.
Prendiamo un primo valore di tentativo $s_1$ e risolviamo l’equazione differenziale con quel valore di tentativo.
Supponiamo che questo primo valore di tentativo ci porti "al di sotto" della nostra soluzione all'istante finale.
Prendiamo un secondo valore di tentativo $s_2$ e supponiamo che questo secondo valore ci porti, invece, "al di sopra" della nostra soluzione all'istante finale.
Ci è capitato un valore sopra ed uno sotto la soluzione in $ t = iTf $, pertanto deduciamo che il nostro "tiro giusto" sarà un valore intermedio tra $s_1 $ ed $s_2 $ , quindi prendiamo un valore $s_3 = (s_1 + s_2) / 2 $
Se in corrispondenza di questo valore $s_3$ la soluzione calcolata all'istante finale è molto prossima allo zero, allora abbiamo concluso; altrimenti applicheremo il metodo di bisezione iterativamente.
Se, ad esempio, con questo valore $s_3$ dovessimo essere ancora al di sopra della soluzione all'istante finale, allora utilizzeremo di nuovo la ricerca con bisezione imponendo che $ s_4 = (s_1 + s_3) / 2 $.
Chiaramente, come facciamo a scegliere inizialmente proprio due valori tali che il primo ci porti "sotto" ed il secondo "sopra" la soluzione all'istante finale?
L'idea è questa; inviamo alla function che implementa lo shooting alcuni parametri $s$ di prova, e andiamo a visualizzare, in corrispondenza di ognuno di questi valori di prova, il valore della soluzione THETA all'istante finale.
In questa figura possiamo apprezzare, in maniera preliminare, il parametro di shooting in corrispondenza del quale la soluzione THETA, all'istante finale, è nulla.
Per inizializzare lo shooting allora basterà prendere due parametri dal grafico, uno che ci faccia andare "sotto" la nostra soluzione ed un altro che ci faccia andare "sopra" la nostra soluzione.
Verifichiamo che questi due valori ci facciano andare "sotto" e "sopra" la soluzione
Se non apparirà un errore, allora i due valori iniziali di shooting saranno esatti e possiamo applicare il metodo di bisezione fintantoché la soluzione in $ t = iTf $ non sia prossima allo zero.
Allora abbiamo in MatLab
Volendo scrivere per bene l'ultimo script si ha
che porta alla seguente soluzione
Supponiamo di essere sotto l'ipotesi di piccoli angoli, in modo tale da approssimare $sin(\theta) = \theta$.
Supponiamo che l'angolo iniziale sia pari a $\theta_0 = 1° $ e che l'angolo finale sia pari a $\theta_f = 0°$.
Supponiamo ancora di avere un pendolo lungo $0.51 m$; allora il suo periodo sarà, sotto l'ipotesi di piccoli spostamenti, $ T = 1.4s$.
Osserviamo poi che siamo in presenza di un problema ai limiti. Abbiamo, infatti, il valore dell'angolo all'istante iniziale e finale. Ci occorre, per risolvere il problema, il valore della derivata prima all'istante iniziale, per cui dobbiamo implementare una procedura di shooting.
Siccome MatLab gestisce solo equazioni di primo grado, dobbiamo trasformare l'equazione del pendolo in un sistema di due equazioni di primo grado.
Ecco lo script
function K = pendolo(theta) global g L theta = theta(:); K(1) = theta(2); K(2) = -(g/L)*(theta(1)); K = K(:)'; end
La funzione riceve in ingresso una matrice theta di due colonne ed una riga, per ogni istante di tempo. Nella prima colonna sono situati tutti i valori della funzione, mentre nella seconda colonna sono conservati tutti i valori della derivata prima della funzione.
Ad esempio, nell'istante di tempo iniziale, la prima colonna conterrà il valore dell'angolo iniziale, mentre la seconda colonna quello della sua derivata prima.
Per gestirlo in MatLab usiamo tuttavia un vettore colonna, ma la sostanza non cambia.
Questa prima function calcola le pendenze, sotto forma di coefficienti K, che servono per il metodo RK di 4 passi.
Creiamo allora lo script che, ricevuto in ingresso il vettore tempo e le condizioni iniziali, implementi il metodo di RK4 utilizzando le K trovate nello script precedente.
Per il metodo di RK4 si può utilizzare il classico array di Butcher
$ \theta_(n+1) = \theta_n + (h)/6*(K1 + 2K2 + 2K3 + K4) $
dove
$ a_(ij) = [ ( 0 , 0, 0, 0),( 0, 1/2, 0, 0),( 0, 0, 1/2, 0),( 0, 0, 0, 1) ] $
$ b = (1/6, 1/3,1/3,1/6) $
$ c = (0, 1/2, 1/2, 1) $
Il metodo di RK4 è implementato da un ciclo for.
La soluzione la mettiamo in una matrice THETA di iTf righe e 2 colonne; come detto sopra, nella prima colonna abbiamo la funzione, mentre nella seconda abbiamo la sua derivata prima (in questo caso spazio e velocità).
Questo script funziona in questo modo;
Il primo istante di tempo è pari alla condizione iniziale; quindi riempio immediatamente la prima riga della matrice THETA con le I.C.
Per gli istanti successivi calcoliamo la funzione THETA, in ogni istante temporale, grazie al metodo RK4.
Supponendo di trovarci all'istante it = 2;
-Prima invio alla function "pendolo" il vettore theta che contiene la condizione iniziale, in modo tale da calcolare la pendenza nel punto iniziale (K1).
-La pendenza nel punto iniziale viene utilizzata per giungere fino al punto $ t + (dt)/2 $;
theta2 = theta + (dt/2)*K1;
In questo punto si calcola poi una nuova pendenza (K2), grazie alla function"pendolo", e la si "riporta" di nuovo al punto iniziale
theta3 = theta + (dt/2)*K2;
Su questa nuova pendenza K2 ci si sposta fino a raggiungere il punto $ t + (dt)/2 $ dove si calcolerà una nuova pendenza (K3), sempre grazie alla function "pendolo", che si "riporterà" al punto iniziale.
Su quest'altra pendenza K3 ci si sposta fino a raggiungere il punto $ t + dt $ dove si calcolerà una nuova pendenza (K4), ancora grazie alla function "pendolo", pendenza che, insieme alle altre tre precedentemente trovate, verrà utilizzata per aggiornare all'istante temporale successivo il valore della THETA.
Veniamo adesso allo shooting.
Per capire questo script dobbiamo fare un piccolo riepilogo.
Noi stiamo risolvendo l'equazione del secondo ordine, tradotta in un sistema di EDO, con un metodo RK4.
Se avessimo tutte e due le condizioni iniziali, problema risolto!
Invece a noi manca la seconda condizione iniziale che va ricercata con il metodo di shooting.
Assegniamo quindi una certa condizione iniziale, $ s $, sperando che sia giusta. Ma come facciamo a sapere che la condizione iniziale "sparata" sia esatta?
Semplicemente, confrontando il valore della funzione in $ t = iTf $, che viene fuori grazie alla risoluzione con
RK4 più seconda condizione iniziale ricavata con shooting ( THETA(end,1) ), e il valore reale che la funzione
assume in $ t = iTf $ .
Adesso, poiché abbiano salvato le soluzioni del sistema di EDO in una matrice dove ogni i-esima riga corrisponde ad un istante temporale i-esimo ed ogni colonna j-esima rappresenta la rispettiva derivata di ordine j-1, allora diciamo "prendi il valore, contenuto nella matrice delle soluzioni THETA, alla riga ultima e pima colonna e confrontalo col valore che so che la funzione assume all'istante finale".
In questo caso, il valore all'istante finale è nullo ($\theta_f = 0° $), quindi dobbiamo ricercare quel parametro di shooting che rende nullo il valore della funzione THETA in $ t = iTf $.
Per lo shooting useremo un metodo di bisezione.
Prendiamo un primo valore di tentativo $s_1$ e risolviamo l’equazione differenziale con quel valore di tentativo.
Supponiamo che questo primo valore di tentativo ci porti "al di sotto" della nostra soluzione all'istante finale.
Prendiamo un secondo valore di tentativo $s_2$ e supponiamo che questo secondo valore ci porti, invece, "al di sopra" della nostra soluzione all'istante finale.
Ci è capitato un valore sopra ed uno sotto la soluzione in $ t = iTf $, pertanto deduciamo che il nostro "tiro giusto" sarà un valore intermedio tra $s_1 $ ed $s_2 $ , quindi prendiamo un valore $s_3 = (s_1 + s_2) / 2 $
Se in corrispondenza di questo valore $s_3$ la soluzione calcolata all'istante finale è molto prossima allo zero, allora abbiamo concluso; altrimenti applicheremo il metodo di bisezione iterativamente.
Se, ad esempio, con questo valore $s_3$ dovessimo essere ancora al di sopra della soluzione all'istante finale, allora utilizzeremo di nuovo la ricerca con bisezione imponendo che $ s_4 = (s_1 + s_3) / 2 $.
Chiaramente, come facciamo a scegliere inizialmente proprio due valori tali che il primo ci porti "sotto" ed il secondo "sopra" la soluzione all'istante finale?
L'idea è questa; inviamo alla function che implementa lo shooting alcuni parametri $s$ di prova, e andiamo a visualizzare, in corrispondenza di ognuno di questi valori di prova, il valore della soluzione THETA all'istante finale.
clc; clear all; close all;
global g L
g = 9.81; L = 0.51;
Tf = 1.4;
iTf = 150;
t = linspace(0, Tf, iTf);
svect = linspace(-500,500,100);
for is=1:length(svect)
s = svect(is);
f(is) = shooting(s);
end
figure(1);
plot (svect, f, 'r.-'); grid on;
title('Zeri della funzione f = THETA(end,1)');
xlabel('Parametro di shooting');
ylabel('f');
In questa figura possiamo apprezzare, in maniera preliminare, il parametro di shooting in corrispondenza del quale la soluzione THETA, all'istante finale, è nulla.
Per inizializzare lo shooting allora basterà prendere due parametri dal grafico, uno che ci faccia andare "sotto" la nostra soluzione ed un altro che ci faccia andare "sopra" la nostra soluzione.
s1 = -100; s2 = 400;
Verifichiamo che questi due valori ci facciano andare "sotto" e "sopra" la soluzione
ys1 = shooting(s1);
ys2 = shooting(s2);
if ys1*ys2 > 0
error('Ridefinisci i parametri iniziali di shooting!');
else
disp('Ok!')
end
Se non apparirà un errore, allora i due valori iniziali di shooting saranno esatti e possiamo applicare il metodo di bisezione fintantoché la soluzione in $ t = iTf $ non sia prossima allo zero.
Allora abbiamo in MatLab
ym = 1;
Tol = 1e-10;
i = 0;
while abs(ym) > Tol
i = i+1;
xm = (s1+s2)/2;
ym = shooting(xm);
if ys1*ym < 0
s2 = xm;
else
s1 = xm;
end
endVolendo scrivere per bene l'ultimo script si ha
che porta alla seguente soluzione
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo