Eulero implicito

[sophii]

Ho implementato il codice relativo a questo problema:
Si risolva il sistema di ODEs

$\{ A'(t)=-2*\alpha(t)*A(t), \alpha'(t)=A(t)^2+\omega(t)^2-\alpha(t)^2-1, \omega'(t)=-2+(\alpha(t)+A(t)*\omega(t)):}$

con dato iniziale

$\{ A(0)=0.5, \alpha(0)=2, \omega(0)=10:}$

con il metodo di Eulero Implicito fino ad un tempo finale $t* = 15$, producendo un grafico della quantità
$E(t) = (A(t)^2+\alpha(t)^2+\omega(t)^2+1)/(2*A(t))$. Si confrontino le soluzioni ottenute usando 300 o 900 timesteps.
% EULERO IMPLICITO, THETA = 1

Codice:

clc
clear all
close all

mrange = [300 , 900];
tstar = 15;
y0 = [0.5; 2; 10

f = @(y,yn) [-2*y(2)*y(1); y(1)^2+y(3)^2-y(2)^2-1; -2*(y(2)+y(1))*y(3)];
jf = @(y,yn) [-2*y(2) , -2*y(1) , 0; 2*y(1) , -2*y(2) , 2*y(3); -2*y(3) , -2*y(3) , -2*(y(2)+y(1))];
F = @(yn,y,k) y-yn-k*f(y);
JF = @(y,k) eye(3)-k*jf(y);

% Soluzione con 300 timesteps
m = mrange(1);
k = tstar/m;
y = y0;
tol = k^2/100;
for j = 1:m
    % Risolvo con Newton poichè è implicito
    y(:,j+1) = y(:,j);
    delta = -JF(y(:,j+1),k)\F(y(:,j),y(:,j+1),k);
    % calcolo della nuova y_(n+1)
    while norm(delta,Inf)>tol
        y(:,j+1) = y(:,j+1)+delta;
        delta = -JF(y(:,j+1),k)\F(y(:,j),y(:,j+1),k);
    end
end

E_300 = (y(1,:).^2+y(2,:).^2+y(3,:).^2+1)./(2*y(1,:)); % si tratta della E del testo!
t_300 = linspace(0,tstar,m+1);
y_300 = y;

% Soluzione con 900 timesteps
m = mrange(2);
k = tstar/m;
y = y0;
tol = k^2/100;
for j = 1:m
    y(:,j+1) = y(:,j);
    delta = -JF(y(:,j+1),k)\F(y(:,j),y(:,j+1),k);
    while norm(delta)>tol
        y(:,j+1) = y(:,j+1)+delta;
        delta = -JF(y(:,j+1),k)\F(y(:,j),y(:,j+1),k);
    end
end
E_900 = (y(1,:).^2+y(2,:).^2+y(3,:).^2+1)./(2*y(1,:));
t_900 = linspace(0,tstar,m+1);
y_900 = y;

% Plottiamo le tre componenti separatamente con 300 e 900 timesteps
subplot(1,3,1)
plot(t_300,y_300(1,:),'x',t_900,y_900(1,:),'o')
subplot(1,3,2)
plot(t_300,y_300(2,:),'x',t_900,y_900(2,:),'o')
subplot(1,3,3)
plot(t_300,y_300(3,:),'x',t_900,y_900(3,:),'o')

% Plottiamo l'energia
figure(2)
plot(t_300,E_300,'x')
figure(3)
plot(t_900,E_900,'x')

% Plottiamo tutte le componenti insieme
figure(4)
plot(t_300,y_300,'x', t_900,y_900,'o')



Mi da errore nella riga 20,ma non capisco perchè
Grazie in anticipo per l'eventuale aiuto.

[feddy]

N.B:

Questa risposta fa riferimento al post che hai modificato circa il confronto tra equazione del pendolo, e equazione del pendolo linearizzato.

L'impostazione del codice va bene. La linearizzazione alla fine consiste in uno sviluppo di Taylor arrestato al prim'ordine del seno, e ti permette di evitare il metodo di Newton.

Il tuo errore è nell'utilizzo del metodo di Newton, hai scambiato i ruoli di $y_{n+1}, y_{n}$ nella funzione $F$ dentro a Newton, che presumo sia quella che vuoi azzerare E' per questo motivo che ti viene il warning "Matrix singular".

Il tuo codice corretto è questo:

Codice:

clear all
close all

m = 100;
l = 1;
g = 9.81;
%accelerazione di gravità)
tstar = pi*sqrt(l/g);
k = tstar/m;
theta0 = pi/4;
y = [theta0; 0];
% inizializzo le funzioni per pendolo non linearizzato
f = @(y) [y(2); -(g/l)*sin(y(1))];
jf = @(y) [0 1 ; 0 -(g/l)*cos(y(1))];
% METODO DEI TRAPEZI
F = @(y,yn) y-yn-k/2*(f(yn)+f(y));
JF = @(y) eye(length(y)) - k/2*jf(y);


% inizializzo le funzioni per il pendolo linearizzato
A = [0 1 ; -g/l, 0];
y2 = y(:,1);
I = eye(length(y2));

tol = k^2/100;
t = 0;

for i = 1:m
    % risoluzione non linearizzato: applico Newton
    t =  t+k;
    % guess iniziale
    y(:,i+1) = y(:,i);
   
    delta = -JF(y(:,i+1))\F(y(:,i+1),y(:,i));
    while norm(delta)>tol
        y(:,i+1) = y(:,i+1)+delta;
        delta = -JF(y(:,i+1))\F(y(:,i+1),y(:,i));
    end
    % risoluzione linearizzato
    y2(:,i+1) = (I-k/2*A)\((I+k/2*A)*y2(:,i));
end       
       
% Trovare il numero minimo di passi affichè la soluzione attraverso Eulero
% disti da theta(tstar) meno di 10^-1

% Guess iniziali
tol = 1;
meul = 10;

while tol > 1e-2
    keul = tstar/meul;
    yeul = [theta0;0];
    for j = 1:meul
        % Applico Eulero esplicito
        yeul(:,j+1) = yeul(:,j)+keul*A*yeul(:,j);
    end
    % Trovo nuova tolleranza e aggiorno meul
    tol = abs(yeul(1,end)-y2(1,end));
    meul = meul+1;
 end
meul = meul-1;
   
subplot(1,2,1)
plot(0:k:tstar,y(1,:),'x',0:k:tstar,y2(1,:),'o')
subplot(1,2,2)
plot(0:k:tstar,y2(1,:),'o',0:keul:tstar,yeul(1,:),'.')


[feddy]

L'errore è di sintassi... non hai definito correttamente $y_0$, ti manca una quadra alla fine e per il resto va, anche se i plot sono un po' incasinati. Quello che importa credo sia mostrare come un maggior numero di time-steps produca una soluzione qualitativamente migliore.

Per il resto il codice mi pare corretto. Nella $F$ da azzerare con Newton, che tu definisci

Codice:

F = @(yn,y,k) y-yn-k*f(y);

metti una dipendenza da $k$ che in realtà non ti serve visto che il passo è costante e $k$ viene definito una volta per tutte come $k=\frac{t_f - t_i}{m}$. Appesantisci solo la tua funzione.

[sophii]

Capito.. Grazie mille

[feddy]

prego!