metodo di Shooting e Runge-Kutta

[sophii]

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<x<3$.
Risolvere con il metodo di shooting sapendo che la soluzione analitica è: $u(x)=x^2+16/x$.
Se possibile risolvere anche con il metodo di runge-kutta.

[feddy]

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

[sophii]

Ho provato e questo è quello che ho ottenuto:

Codice:


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

[feddy]

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

[feddy]

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


Scelta del metodo
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.

Bisezione

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

Codice:

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

[feddy]

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è

Codice:

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

Codice:

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:

Codice:

u0=linspace(a,b,m)';

Una sua possibile implementazione, facendo riferimento al tuo codice, è questa:

Codice:

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.

Codice:

  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


[sophii]

Tutto chiaro..Grazie

[feddy]

prego