% ----------------------------------------
% Hauptfile zu
% dglsyst_MND2(...)
% Methode der Taylorreihe
% im Vergleich zur exakten Lösung
% Serie7, Aufgabe 6
% n = 3 A = 3 mal 3 Matrix
% März 2014
% --------------------------------------------------

% AB:
alpha = 1.0;
beta = 0.4;
gamma = 2;
x0 = [alpha; beta; gamma]

t0 = 0.0;
tend = 8.0;

% ng = Anzahl Schritte mit der Methode der Taylorreihe
ng = 10*tend;
delta_t = (tend - t0)/ng;
t = t0:delta_t:tend;

x = [];
x = [x x0];

p = input('Anzahl Terme der Taylorreihe p = ')

% Methode der Taylorreihe
% Entwicklung bis und mit Term p-ter Ordnung

for k = 2:length(t)
    c = zeros(length(x0), p);
    c(:,1) = feval('dglsyst_2014', t(k-1), x(:,k-1));
    for j = 2:p
     c(:,j) = 1/j*feval('dglsyst_2014', t(k-1), c(:,j-1));
    end
    xk = x(:,k-1);
    for j= 1:p
        xk = xk + delta_t^j*c(:,j);
    end
     x = [x, xk];
    
end
    
figure(1)
hold off
plot(t, x(1,:), 'r.-')
hold all
plot(t, x(2,:), '.-')
hold all
plot(t, x(3,:), '.-')
grid
legend('x1(t) = ', 'x2(t) = ', 'x3(t) = ')
title('Komponenten einzeln')
xlabel(['Schrittweite h = ', num2str(delta_t)])

pause

figure(2)
hold off
plot(x(1,:), x(2,:), 'b.-')
grid
    
title('Phasenporträt x1 x2')

pause

figure(3)
hold off
plot(x(2,:), x(3,:), 'g.-')
grid
title('Phasenporträt x2 x3')
pause

figure(4)
hold off
plot3(x(1,:), x(2,:),x(3,:), 'k.-')
grid

% Vergleich mit der exakten Lösung

% A = [-1 1;-3 -3];
% [T,D] = eig(A);
% 
% alpha1 = real(D(1))
% beta1 = imag(D(1))
% 
% u = real(T(:,1));
% w = imag(T(:,1));
% 
% % Bestimmung von a1 und b1:
% 
% B = [u -w];
% 
% c = B\x0;        % a1 = c(1), b1 = c(2)
% a1 = c(1);
% b1 = c(2);
% 
% % Bestimmung der LK
% 
% n = length(t);
% xe = [];
% 
% for k=1:n
%     xk = (a1*cos(beta1*t(k)) - b1*sin(beta1*t(k)))*u ...
%                    -(a1*sin(beta1*t(k)) + b1*cos(beta1*t(k)))*w;
%     xk = exp(alpha1*t(k))*xk;
%     xe = [xe xk];
% end
% 
% figure(3)
% hold off
% plot(x(1,:), x(2,:), 'b.-')
% hold on
% plot(xe(1,:), xe(2,:), 'k.-')
% grid
% legend('Taylorreihe', 'exakt', 2)
% xlabel(['Anzahl Terme in der Entwicklung: p = ', num2str(p)])
% title('Phasenporträt')
% 
% pause
% 
% % Differenzen:
% disp('Differenzen: x - xe')
% x - xe
% % absoluter Fehler:
% absF = zeros(1, length(t));
% relF = zeros(1, length(t));
% 
% for k = 1:length(t)
%  absF(k) = norm(x(:,k) - xe(:,k));
%  relF(k) = absF(k)/norm(xe(:,k));
% end
% 
% max(absF)
% 
% figure(4)
% hold off
% semilogy(t, absF, 'm.-')
% hold on
% %semilogy(t, relF, 'k.-')
% grid
% title('Fehler')
% %legend('absolute Fehler', 'relative Fehler', 2)
% legend('absolute Fehler')
% xlabel(['Schrittweite h = ', num2str(delta_t)])
% 
