% Num. Math. - FS 2013 - Serie 5 Aufgabe 2 - Implicit midpoint rule / Heun


function Implicitmidpoint
clear all;close all;
%Anfangsbedingung
y0=[0.5;0;0;sqrt(3)];
%Endzeitpunkt
tf=500;
%Vektor mit versch. Schrittweiten fuer beide Verfahren
h_vec=[1/5; 1/10;1/100];
%Iteration ueber die Schrittweiten
for i=1:length(h_vec)
 h=h_vec(i); %Schrittweite
 
 N=...% !!!!!!!    %Anzahl Schritte
 
 y_heun=zeros(length(y0),N+1); %Initializierung der Loesungs'matrix' 
                               %fuer Heun
 y_imp_mp=zeros(length(y0),N+1); %Initializierung der Loesungs'matrix'
                                 % fuer impl. Mittelpunktsregel
 y_heun(:,1)=y0; %Initialisierung der Anf.bedingung
 y_imp_mp(:,1)=y0; %Initialisierung der Anf.bedingung
 %Iteration in der Zeit 
 for k=1:N
   y_heun(:,k+1)=heun_schritt(@yprime,(k-1)*h,y_heun(:,k),h);
 
   y_imp_mp(:,k+1)=impl_mp_step(..);% !!!!!
  
 end
 %Plotten der ersten beiden Komponenten 
 subplot(3,2,2*i-1);
 plot(y_imp_mp(1,:),y_imp_mp(2,:));
 axis equal; axis([-3 2 -1.5 2]);
 title(['Implizite Mittelpunktsregel mit h=',num2str(h),...
     ', Bahn von t=0 bis t=',num2str(tf)])  
 subplot(3,2,2*i);
 plot(y_heun(1,:),y_heun(2,:));
 axis equal; axis([-3 2 -1.5 2]);
 title(['Heun mit h=',num2str(h),...
     ', Bahn t=0 bis t=',num2str(tf)]) 
end
%Speichern des Diagramms
print('-depsc','pl5_2.eps')
end

function[y_neu]=heun_schritt(fct,t,y,h)
 f_=fct(t,y);
 y_neu=y+h*f_;
 y_neu=y+h/2*(f_+fct(t+h,y_neu));
end

function[y_neu]=impl_mp_schritt(fct,t,y,h)
%Parameter fuer Fixpunktiteration
tol=1e-6;
maxit=100;

y_neu=...      % !!! Startwert fuer Fixpunktiteration: Euler

y_vorh_schritt=y_neu;
for i=1:maxit
       
    y_neu=... %!!! Fixpunktiteration 
  
    %Abbruchkriterium
    if(norm(y_neu-y_vorh_schritt)<tol+tol*norm(y_neu))
        return
    end
    y_vorh_schritt=y_neu;
end
disp('Fixpunktiteration hat nicht konvergiert')
end

function[f_]=yprime(t,y)
denom=norm(y(1:2))^3;
f_=[y(3);y(4);-y(1)/denom;-y(2)/denom];
end
