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.