Metodo di Shooting e Runge-Kutta
Ciao a tutti. Non riesco più ad andare avanti con il seguente problema , faccio fatica con le condizioni al bordo. Il sistema è:
$ (u''(x)=1/8*(32+2*x^3-u(x)*u'(x)), u(1)=17 , u(3)=43/3): $ con $1
Se possibile risolvere anche con il metodo di runge-kutta.
$ (u''(x)=1/8*(32+2*x^3-u(x)*u'(x)), u(1)=17 , u(3)=43/3): $ con $1
Se possibile risolvere anche con il metodo di runge-kutta.
Risposte
In questo caso metodo di shooting significa trasformare il problema ai limiti in un sistema di ODE del prim'ordine e poi utilizzare un metodo di bisezione, o comunque di ricerca di zeri, in modo opportuno per ricavare il valore delle condizioni iniziale $y_2(0)$.
Cioé il problema ai limiti diventa una famiglia di problemi ai valori iniziali che dipenderà da $r$:
\( \begin{cases} y_1'=y_2 \\ y_2'=\frac{1}{8} (32 + 2t^3 - y_1 \cdot y_2) \\ y_1(1)=17 \\ y_2(1)=r \end{cases} \)
che puoi risolvere col metodo che preferisci, nel tuo caso un Runge Kutta, ma si potrebbero anche usare metodi più semplici.
Una volta trovata la soluzione, dovrai implementare uno schema di tipo bisezione per azzerare la funzione $F(r)=\bar{y} - u_b$. Il significato è trovare un valore di $r$ (seconda componente della condizione iniziale) tale che la funzione valutata nell'estremo (indicato con $\bar(y)$) abbia il valore della condizione al bordo del problema ai limiti iniziale $u_b$. Volendo essere più corretti bisognerebbe scrivere $\bar(y)=\bar(y)(r)$
Cioé il problema ai limiti diventa una famiglia di problemi ai valori iniziali che dipenderà da $r$:
\( \begin{cases} y_1'=y_2 \\ y_2'=\frac{1}{8} (32 + 2t^3 - y_1 \cdot y_2) \\ y_1(1)=17 \\ y_2(1)=r \end{cases} \)
che puoi risolvere col metodo che preferisci, nel tuo caso un Runge Kutta, ma si potrebbero anche usare metodi più semplici.
Una volta trovata la soluzione, dovrai implementare uno schema di tipo bisezione per azzerare la funzione $F(r)=\bar{y} - u_b$. Il significato è trovare un valore di $r$ (seconda componente della condizione iniziale) tale che la funzione valutata nell'estremo (indicato con $\bar(y)$) abbia il valore della condizione al bordo del problema ai limiti iniziale $u_b$. Volendo essere più corretti bisognerebbe scrivere $\bar(y)=\bar(y)(r)$
Ho provato e questo è quello che ho ottenuto:
Però ci sono ancora dei buchi: non so se le correzioni per le condizioni al bordo sono giuste (non penso, in caso non so come fare) e manca il vettore iniziale per applicare newton (penso di poterlo usare in questo caso).
clc
clear all
close all
a=1;
b=3;
mrange = [101,201,251,501,641,751];
counter = 0;
for i=1:length(mrange)
m = mrange(i);
counter = counter+1;
x = linspace(a,b,m)';
h = (b-a)/(m-1);
tol = 1e-6;
A = toeplitz(sparse([1,1],[1,2],[-2,1]/h^2,1,m));
% correzione per i bordi
A(1,1:2) = [1,0]; %???
A(m,m-1:m) = [0,1];%???
%matrice delle derivate prime
B=toeplitz(sparse(1,2,-1/(2*h),1,m),sparse(1,2,1/(2*h),1,m));
B(1,1:2)=[0,0]; %%???
B(m,m-1:m)=[0,0];%???
rhs= %???
F = @(u) A*u+(1/8).*u.*(B*u)-4-(1/4)*x.^3+rhs;
JF =@(u) A+(1/8)*diag(B*u)+(1/8)*diag(u)*B;
% NEWTON
% Vettore iniziale
u0 = % ???
u = u0;
delta = -JF(u)\F(u);
while (norm(delta,inf)>tol)
u=u+delta;
delta=-JF(u)\F(u);
end
uesatta=x.^2+16./x;
% Calcolo l'errore
err(counter)=norm(u-uesatta,inf);
end
figure(1)
loglog(mrange,err,'*',mrange,mrange.^(-2),'r');
title('ANDAMENTO dell' ERRORE')
legend('errore','retta di pendenza -2')
xlabel('Asse x')
ylabel('Asse y')
grid on
Però ci sono ancora dei buchi: non so se le correzioni per le condizioni al bordo sono giuste (non penso, in caso non so come fare) e manca il vettore iniziale per applicare newton (penso di poterlo usare in questo caso).
Un metodo di shooting è un problema ai valori iniziali, quello che hai messo è solo un problema ai limiti.
Sono due cose diverse.
Nel primo caso hai due condizioni iniziali nel primo bordo. Nel secondo due valori agli estremi dell'intervallo.
Nel primo caso usi un metodo come Runge Kutta. Nel secondo in genere differenze finite
Sono due cose diverse.
Nel primo caso hai due condizioni iniziali nel primo bordo. Nel secondo due valori agli estremi dell'intervallo.
Nel primo caso usi un metodo come Runge Kutta. Nel secondo in genere differenze finite
Nel caso di un metodo di shooting:
Definizione della funzione $f$ del problema. In questo caso sarà $f: (t,y) \rarr RR^2$, soddisfacente a tutte le opportune condizioni (esistenza, unicità)
Ovviamente la matrice jacobiana è \( J_f=\begin{matrix} 0 & 1 \\ -\frac{y_2}{8} & -\frac{y_1}{8} \end{matrix} \)
Discretizzare l'intervallo temporale con $n$ passi temporali di ampiezza uguale a $h$, ossia $h=\frac{b-a}{n}$
[size=150]Scelta del metodo[/size]
Per semplicità supponiamo di usare eulero implicito.
Lo schema in questo caso è
Per trovare il valore di $y_{n+1}$ devi risolvere un sistema non lineare, e quindi usare il metodo di Newton.
[size=150]Bisezione[/size]
Per prima cosa serve individuare un intervallo plausibile. L'incognita $r$ altro non è che la pendenza iniziale della derivata, pertanto, senza fare studi qualitativo di alcun genere, prendiamo come intervallo $[-20,10]$.
Applichiamo il metodo di Eulero implicito con dato iniziale $y_{0}=[17;a]$ e otteniamo un certo valore al bordo $x=3$, che chiamiamo $y_a$. Quello che si vuole ottenere è che vogliamo azzerare è $y_a - \frac{43}{3}$, e chiamiamo questa quantità $f_a$.
Analogamente a quanto fatto sopra, applichiamo Eulero implicito con dato iniziale $[17;b]$. Otteniamo un valore a $x=3$ che vale $y_b$, e come prima definiamo $f_b= y_b - \frac{43}{3}$.
Ora devi solo fare bisezione su $a$ e $b$. In un certo senso, rispetto al metodo di bisezione a cui si è abituati, il ruolo dell'intervallo è giocato da $[a,b]$, mentre i valori della funzione da azzerare sono rappresentati da $f_a, f_b$.
Ripeti questo procedimento fino a che la $|b-a| < \text{Tol} \quad && \quad \text{iter}<\text{MaxIter}$. (criterio d'arresto).
Definizione della funzione $f$ del problema. In questo caso sarà $f: (t,y) \rarr RR^2$, soddisfacente a tutte le opportune condizioni (esistenza, unicità)
Ovviamente la matrice jacobiana è \( J_f=\begin{matrix} 0 & 1 \\ -\frac{y_2}{8} & -\frac{y_1}{8} \end{matrix} \)
Discretizzare l'intervallo temporale con $n$ passi temporali di ampiezza uguale a $h$, ossia $h=\frac{b-a}{n}$
[size=150]Scelta del metodo[/size]
Per semplicità supponiamo di usare eulero implicito.
Lo schema in questo caso è
$y_{n+1}=y_{n} + h \cdot f(t_{n+1},y_{n+1})$
.Per trovare il valore di $y_{n+1}$ devi risolvere un sistema non lineare, e quindi usare il metodo di Newton.
[size=150]Bisezione[/size]
Per prima cosa serve individuare un intervallo plausibile. L'incognita $r$ altro non è che la pendenza iniziale della derivata, pertanto, senza fare studi qualitativo di alcun genere, prendiamo come intervallo $[-20,10]$.
Applichiamo il metodo di Eulero implicito con dato iniziale $y_{0}=[17;a]$ e otteniamo un certo valore al bordo $x=3$, che chiamiamo $y_a$. Quello che si vuole ottenere è che vogliamo azzerare è $y_a - \frac{43}{3}$, e chiamiamo questa quantità $f_a$.
Analogamente a quanto fatto sopra, applichiamo Eulero implicito con dato iniziale $[17;b]$. Otteniamo un valore a $x=3$ che vale $y_b$, e come prima definiamo $f_b= y_b - \frac{43}{3}$.
Ora devi solo fare bisezione su $a$ e $b$. In un certo senso, rispetto al metodo di bisezione a cui si è abituati, il ruolo dell'intervallo è giocato da $[a,b]$, mentre i valori della funzione da azzerare sono rappresentati da $f_a, f_b$.
xm=(a+b)/2; %punto medio
y(:,1)=[17,xm];
%[Eulero implicito]
fm=y(1,end)-43/3;
if sign(fa)*sign(fm)<=0 %bisezione
b=xm;
fb=fm;
else
a=xm;
fa=fm;
end
end
Ripeti questo procedimento fino a che la $|b-a| < \text{Tol} \quad && \quad \text{iter}<\text{MaxIter}$. (criterio d'arresto).
Nel caso di un problema ai limiti:
\( \bullet \) Per le condizioni ai bordi ti consiglio di scrivere su un foglio la prima e l'ultima riga, in modo da avere chiaro cosa dover "azzerare", o comunque cosa dover modificare.
Per la condizione a $x=1$ basterebbe azzerare la prima riga della matrice che discretizza la derivata prima e modificare la prima riga della matrice che discretizza la derivata seconda lasciando solo l'entrata (1,1) pari a $1$.
Cioè
A questo punto, definiamo il termine noto $b= -4 - \frac{x^3}{3}$, con la prima e ultima entrata rispettivamente
\( \bullet \) Definiamo la funzione da azzerare $F(u)= Au + \frac{1}{8} (\diag(u) \dot ( B \dot u ) + \vec{b})$
Evidentemente è un funzionale non lineare e dunque come hai visto bisogna usare Newton.
\( \bullet \) Lo parte più difficile per fare lo jacobiano è quella riguardante il termine $\diag(u) \cdot (B \cdot u)$.
Ci sono più modi: il più semplice consiste nello scrivere l'i-esima riga generale di questo prodotto:
Derivando rispetto alla i-1-esima: $-\frac{u_i}{2h}$
Derivando rispetto alle i-esima: $\frac{u_{i+1}-u_{i-1}}{2h}$
Derivando rispetto alla i+1-esima: $\frac{u_i}{2h}$
Quello che risulta è una matrice tridiagonale, che è equivalente a quella scritta da te.
La matrice jacobiana è dunque $J=A + \frac{1}{8} \cdot (diag(B \cdot u) + diag(u) \cdot B)$
\( \bullet \) Il metodo di Newton, essendo iterativo, prevede un guess iniziale. Sotto opportune ipotesi Newton converge quadraticamente. Un buon guess iniziale può essere una retta congiungente i due estremi:
Una sua possibile implementazione, facendo riferimento al tuo codice, è questa:
Nel codice ho mostrato anche che l'ordine di convergenza è effettivamente $2$, mentre per risolvere il sistema non lineare ho utilizzato un criterio d'arresto misto con Newton inesatto (cioè usando sempre la stessa matrice jacobiana). Provando a cambiare il guess iniziale in quello da me suggerito, dovresti vedere che il metodo impiega meno iterazioni (passi da 28 iterate a 25).
Inoltre, usando il metodo di Newton (non inesatto quindi), hai che in 4 iterazioni hai raggiunto la soluzione.
\( \bullet \) Per le condizioni ai bordi ti consiglio di scrivere su un foglio la prima e l'ultima riga, in modo da avere chiaro cosa dover "azzerare", o comunque cosa dover modificare.
Per la condizione a $x=1$ basterebbe azzerare la prima riga della matrice che discretizza la derivata prima e modificare la prima riga della matrice che discretizza la derivata seconda lasciando solo l'entrata (1,1) pari a $1$.
Cioè
A(1,1:2)=[1,0]; B(1,2)=0; A(m,m-1:m)=[0,1]; B(m,m-1)=0;
A questo punto, definiamo il termine noto $b= -4 - \frac{x^3}{3}$, con la prima e ultima entrata rispettivamente
b(1)=- 17/(h^2); b(m)= - 43/(3*h^2);
\( \bullet \) Definiamo la funzione da azzerare $F(u)= Au + \frac{1}{8} (\diag(u) \dot ( B \dot u ) + \vec{b})$
Evidentemente è un funzionale non lineare e dunque come hai visto bisogna usare Newton.
\( \bullet \) Lo parte più difficile per fare lo jacobiano è quella riguardante il termine $\diag(u) \cdot (B \cdot u)$.
Ci sono più modi: il più semplice consiste nello scrivere l'i-esima riga generale di questo prodotto:
$u_i \cdot \frac{u_{i+1} - u_{i-1}}{2h}$
Derivando rispetto alla i-1-esima: $-\frac{u_i}{2h}$
Derivando rispetto alle i-esima: $\frac{u_{i+1}-u_{i-1}}{2h}$
Derivando rispetto alla i+1-esima: $\frac{u_i}{2h}$
Quello che risulta è una matrice tridiagonale, che è equivalente a quella scritta da te.
La matrice jacobiana è dunque $J=A + \frac{1}{8} \cdot (diag(B \cdot u) + diag(u) \cdot B)$
\( \bullet \) Il metodo di Newton, essendo iterativo, prevede un guess iniziale. Sotto opportune ipotesi Newton converge quadraticamente. Un buon guess iniziale può essere una retta congiungente i due estremi:
u0=linspace(a,b,m)';
Una sua possibile implementazione, facendo riferimento al tuo codice, è questa:
u0 = %dato iniziale; AbsTol=1e-9; RelTol=1e-6; delta=-JF(u0)\F(u0); while (norm(delta,inf)> AbsTol + norm(res,inf)*RelTol) u0=u0+delta; delta=JF(u0)\F(u0); end
Nel codice ho mostrato anche che l'ordine di convergenza è effettivamente $2$, mentre per risolvere il sistema non lineare ho utilizzato un criterio d'arresto misto con Newton inesatto (cioè usando sempre la stessa matrice jacobiana). Provando a cambiare il guess iniziale in quello da me suggerito, dovresti vedere che il metodo impiega meno iterazioni (passi da 28 iterate a 25).
Inoltre, usando il metodo di Newton (non inesatto quindi), hai che in 4 iterazioni hai raggiunto la soluzione.
clear all
close all
count=0;
mrange=2.^(4:10)+1;
for m=mrange
count=count+1;
x=linspace (1,3,m)'; % intervallo del problema differenziale
h=(2)/(m-1);%lunghezza dell'intervallo
A = toeplitz(sparse([1,1],[1,2],[-2,1]/(h^2),1,m));
B = toeplitz(sparse(1,2,-1/(2*h),1,m), sparse(1,2,1/(2*h),1,m));
b=zeros(m,1);
%C. AL BORDO
A(1,1:2)=[1/(h^2),0];
B(1,1:2)=[0,0];
A(m,m-1:m)=[0,1/(h^2)];
B(m,m-1:m)=[0,0];
b= -4*ones(m,1) - (1/4)*(x.^3);
b(1,1)=- 17/(h^2);
b(m,1)= - 43/(3*h^2);
F=@(u) A*u + (1/8)*(diag(u)*(B*u))+ b;
JF=@(u) A + 1/8*(diag(B*u) + diag(u)*B);
%PARAMETRI INIZIALI
% u0=ones(m,1);
u0=linspace(1,3,m)';
u=u0;
delta=-JF(u)\F(u);
i=0;
RelTol=1e-6;
AbsTol=1e-9;
while(norm(delta,inf)>RelTol*norm(delta,inf) + AbsTol)
i=i+1;
u=u+delta;
delta=-JF(u0) \ F(u); %newton inesatto, convergenza lineare
end
u=u+delta;
s=@(x) x.^2 + 16./x;
err(count)=norm(u- s(x),inf)
end
figure
plot(x,u,'*',x,x.^2 + 16./x,'r')
title('soluzioni')
legend('sol numerica','sol analitica')
xlabel('x')
ylabel('u(x)')
i %iterate
figure(2)
loglog(mrange,err,'sk',mrange,(err(end))*(mrange/mrange(end)).^-2,'r',mrange,mrange.^(-1),'b')
title('errore')
legend('errore','O(m^-2)','O(m^-1)')
xlabel('m')
ylabel('err')
grid on
Tutto chiaro..Grazie 
prego
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo