Problema metodo shooting.

[Niernen]

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!!

[goffredoTerzo]

Guarda, anche io ho lo stesso problema, cerchiamo di risolverlo insieme entro il 16? Altrimenti non lo diamo calcolo..

Con affetto,

Iustin

[astrolabio95]

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

Codice:

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) $

Testo nascosto, fai click qui per vederlo

Codice:

function THETA = RK4(t,theta0)

iTf = length(t);
dt  = t(2) - t(1);

THETA      = nan(iTf,2);
THETA(1,:) = theta0;

theta = theta0;

for it = 2:iTf
   
    theta1 = theta;
    K1     = pendolo(theta1);
   
    theta2 = theta + (dt/2)*K1;
    K2     = pendolo(theta2);
   
    theta3 = theta + (dt/2)*K2;
    K3     = pendolo(theta3);
   
    theta4 = theta + dt*K3;
    K4     = pendolo(theta4);
   
    theta  = theta + (dt/6)*(K1 + 2*K2 +2*K3 +K4);
   
    THETA(it,:) = theta;
   
end

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

Codice:

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

Codice:

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.

Testo nascosto, fai click qui per vederlo

Codice:

function f = shooting(s)

global g L

g = 9.81; L = 0.51;

Tf  = 1.4;
iTf = 150;
t   = linspace(0, Tf, iTf);

theta0 = [1 s];

THETA = RK4(t,theta0);

f = THETA(end,1);

end


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.

Codice:

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.

Codice:

s1 = -100;
s2 = 400;

Verifichiamo che questi due valori ci facciano andare "sotto" e "sopra" la soluzione

Codice:

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

Codice:

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
end

Volendo scrivere per bene l'ultimo script si ha

Testo nascosto, fai click qui per vederlo

Codice:

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');

s1 = -100;
s2 = 400;

ys1 = shooting(s1);
ys2 = shooting(s2);

if ys1*ys2 > 0
    error('Ridefinisci i parametri iniziali di shooting!');
else
    disp('Ok!')
end

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
end
disp(['Il parametro di shooting è ', num2str(xm), ' trovato dopo ',...
    num2str(i), ' iterazioni.'])

theta0 = [1 xm];

THETA = RK4(t,theta0);

figure(2);
plot(t, THETA(1:iTf,1), 'b'); grid on;
title('Equazione del pendolo semplice. \theta_0 = 1°');
xlabel('t = sec');
ylabel('\theta = deg');


che porta alla seguente soluzione