function A29

%harmonischer Oszillator
omega=2*pi;
A=[0,1;-omega^2,0];
x0=[1,0]';
T=5;
N=1000;
[t,y1]=Eexp(A,x0,T,N);
[t,y2]=Eimp(A,x0,T,N);
y=[cos(omega*t);omega*sin(omega*t)];
plot(y1(1,:),y1(2,:),'b',y2(1,:),y2(2,:),'g',y(1,:),y(2,:),'r');
axis square;
legend('Euler explizit','Euler implizit','exakt');

figure;

N=50;
[t,y3]=Mexp(A,x0,T,N);
[t,y4]=Mimp(A,x0,T,N);
[t,y5]=RK4(A,x0,T,N); 
[t,y6]=Esymp(-omega^2,x0,T,N); 
plot(y3(1,:),y3(2,:),'b-*',y4(1,:),y4(2,:),'go',y5(1,:),y5(2,:),'kd',y6(1,:),y6(2,:),'cs-',y(1,:),y(2,:),'r');
axis square;
legend('Mittelpunkt explizit','Mittelpunkt implizit','Runge Kutta','Euler symplektisch','exakt');

T=2;
N=100;
[t,y1]=Eexp(A,x0,T,N);
[t,y2]=Eimp(A,x0,T,N);
[t,y3]=Mexp(A,x0,T,N);
[t,y4]=Mimp(A,x0,T,N);
[t,y5]=RK4(A,x0,T,N);
[t,y6]=Esymp(-omega^2,x0,T,N);
y=[cos(omega*t);omega*sin(omega*t)];


figure;

%Energie
en1=omega^2*y1(1,:).^2+y1(2,:).^2;
en2=omega^2*y2(1,:).^2+y2(2,:).^2;
en3=omega^2*y3(1,:).^2+y3(2,:).^2;
en4=omega^2*y4(1,:).^2+y4(2,:).^2;
en5=omega^2*y5(1,:).^2+y5(2,:).^2;
en6=omega^2*y6(1,:).^2+y6(2,:).^2;
en=omega^2*y(1,:).^2+y(2,:).^2;
plot(t,en1,'b',t,en2,'g',t,en3,'c',t,en4,'y',t,en5,'k',t,en6,'m',t,en,'r');
legend('Euler explizit', 'Euler implizit', 'Mittelpunkt explizit', 'Mittelpunkt implizit',' Runge Kutta','Euler symplektisch','exakt');
title('Gesamtenergie');

figure;
 
%Konvergenzanalyse
T=5;
N=round(logspace(1,5,10));
%Zählvariable
c=1;
%Initialisierung
e1=zeros(1,length(N));
e2=zeros(1,length(N));
e3=zeros(1,length(N));
e4=zeros(1,length(N));
e5=zeros(1,length(N));
for n=N 
%approx. Lösungen
[t,y1]=Eexp(A,x0,T,n);
[t,y2]=Eimp(A,x0,T,n);
[t,y3]=Mexp(A,x0,T,n);
[t,y4]=Mimp(A,x0,T,n);
[t,y5]=RK4(A,x0,T,n);
[t,y6]=Esymp(-omega^2,x0,T,n);
y=[cos(omega*t);omega*sin(omega*t)];
%Energien
en1=omega^2*y1(1,:).^2+y1(2,:).^2;
en2=omega^2*y2(1,:).^2+y2(2,:).^2;
en3=omega^2*y3(1,:).^2+y3(2,:).^2;
en4=omega^2*y4(1,:).^2+y4(2,:).^2;
en5=omega^2*y5(1,:).^2+y5(2,:).^2;
en6=omega^2*y6(1,:).^2+y6(2,:).^2;
exa=omega^2*y(1,:).^2+y(2,:).^2;
e1(c)=max(abs(exa-en1));
e2(c)=max(abs(exa-en2));
e3(c)=max(abs(exa-en3));
e4(c)=max(abs(exa-en4));
e5(c)=max(abs(exa-en5));
e6(c)=max(abs(exa-en6));
c=c+1;
end;
loglog(N,e1,'b',N,e2,'g',N,e3,'c',N,e4,'y',N,e5,'k',N,e6,'m');
legend('Euler explizit', 'Euler implizit', 'Mittelpunkt explizit', 'Mittelpunkt implizit',' Runge Kutta','Euler symplektisch');

  
function [t,y]=Eexp(A,x0,T,N)
% löst AWP für lineare gewöhnliche DGL mit Matrix A, Startwert x0
% bis zur Zeit T mit N Schritten

% Raumdimension
d=length(x0);
% Lösungsvektor vorbereiten
y=zeros(d,N);
% Zeitpunkte
t=linspace(0,T,N);
% Zeitschritt
dt=t(2)-t(1);
% Anfangswert eintragen
y(:,1)=x0;
% Iterationsmatrix
B=eye(d)+dt*A;
for i=2:N
 y(:,i)=B*y(:,i-1);
end;

function [t,y]=Eimp(A,x0,T,N)
% löst AWP für lineare gewöhnliche DGL mit Matrix A, Startwert x0
% bis zur Zeit T mit N Schritten
 
% Raumdimension
d=length(x0);
% Lösungsvektor vorbereiten
y=zeros(d,N);
% Zeitpunkte
t=linspace(0,T,N);
% Zeitschritt
dt=t(2)-t(1);
% Anfangswert eintragen
y(:,1)=x0;
% Iterationsmatrix
B=(eye(d)-dt*A)^(-1);
for i=2:N
 y(:,i)=B*y(:,i-1);
end;   


function [t,y]=Mexp(A,x0,T,N)
% löst AWP für lineare gewöhnliche DGL mit Matrix A, Startwert x0
% bis zur Zeit T mit N Schritten
 
% Raumdimension
d=length(x0);
% Lösungsvektor vorbereiten
y=zeros(d,N);
% Zeitpunkte
t=linspace(0,T,N);
% Zeitschritt
dt=t(2)-t(1);
% Anfangswert eintragen
y(:,1)=x0;
y(:,2)=(eye(d)+dt*A)*x0;
% Iterationsmatrix
B=2*dt*A;
for i=3:N
 y(:,i)=y(:,i-2)+B*y(:,i-1);
end;

function [t,y]=Mimp(A,x0,T,N)
% löst AWP für lineare gewöhnliche DGL mit Matrix A, Startwert x0
% bis zur Zeit T mit N Schritten
 
% Raumdimension
d=length(x0);
% Lösungsvektor vorbereiten
y=zeros(d,N);
% Zeitpunkte
t=linspace(0,T,N);
% Zeitschritt
dt=t(2)-t(1);
% Anfangswert eintragen
y(:,1)=x0;
% Iterationsmatrix
B=(eye(d)-dt/2*A)^(-1)*(eye(d)+dt/2*A);
for i=2:N
 y(:,i)=B*y(:,i-1);
end;   

function [t,y]=RK4(A,x0,T,N)
% löst AWP für lineare gewöhnliche DGL mit Matrix A, Startwert x0
% bis zur Zeit T mit N Schritten
 
% Raumdimension
d=length(x0);
% Lösungsvektor vorbereiten
y=zeros(d,N);
% Zeitpunkte
t=linspace(0,T,N);
% Zeitschritt
dt=t(2)-t(1);
% Anfangswert eintragen
y(:,1)=x0; 
for i=2:N
 x=y(:,i-1);
 k1=A*x;
 k2=A*(x+dt*k1/2);
 k3=A*(x+dt*k2/2);
 k4=A*(x+dt*k3);
 y(:,i)=x+dt*(k1+2*(k2+k3)+k4)/6;
end;

function [t,y]=Esymp(A,x0,T,N)

% Raumdimension
d=length(x0);
m=d/2;
% Lösungsvektor vorbereiten
y=zeros(d,N);
% Zeitpunkte
t=linspace(0,T,N);
% Zeitschritt
dt=t(2)-t(1);
% Anfangswert eintragen
y(:,1)=x0;
B=dt*A;
for i=2:N 
  x=y([1:m],i-1);
  v=y([m+1:end],i-1);
  vn=v+B*x;
  xn=x+dt*vn;
  y(:,i)=[xn;vn];
end;

