% Num. Math. - FS 2013 - Serie 6 Aufgabe 1 - Van der Pol oscillator

 function s6_1 
 close all;
 clear all;
 epsilon =1e-3;
 T=5;
 y0=[1;1];
 t0=0;
 f=@(t,y) %%%
 Df=@(t,y)[0 1;-1/epsilon 1/epsilon*(1-y(2)^2)];
 
 tol = 1e-2;
 h   = 1e-3;
 
 %=============================================== part a
 
 [t,y]=simple_heun(f,t0,y0,h,T);
 
 subplot(2,3,1)
 plot(y(1,:),y(2,:))
 title('Heun')
 axis equal
 axis([-2 2 -3 3])

 %=============================================== part b
 
 [t,y]=richardson_heun(f,t0,y0,h,T,tol);
 
 subplot(2,3,2)
 plot(y(1,:),y(2,:))
 title('Heun (Richardson)')
 axis equal
 axis([-2 2 -3 3])

 subplot(2,3,3)
 semilogy(diff(t))
 title('Schrittweiten (Richardson)')
 n_heun_steps=length(t)-1
 
 %=============================================== part c
 
 [t,y]=trapezmethode(f,t0,y0,h,T,Df);
 
 subplot(2,3,4);
 plot(y(1,:),y(2,:));
 axis equal
 axis([-2 2 -3 3]);
 title('Trapezmethode')
 
 %=============================================== part d
 
 %MATLAB ODE-Solver ODE23S
 %%%
 [t,y]=%%% ode23s benutzen
 subplot(2,3,5)
 plot(y(:,1),y(:,2));
 axis equal
 axis([-2 2 -3 3])
 title('ODE23S')
 
 subplot(2,3,6);
 semilogy(diff(t))
 title('Schrittweiten ODE23S') 
 
 n_ode23s_steps=length(t)-1
 
 %Diagramm speichern
 print('-depsc','pl6_1.eps')     
 end
 
 
 
 
 function[t,y]=simple_heun(f,t0,y0,h,T)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Verfahren von Heun zum Loesen der ODE y'=f(t,y), y(t0)=y0
% f: Funktionshandle auf f(t,y)
% y0: Anfangswerte
% t0: Anfangszeit
% h: Schrittweite
% T: Endzeitpunkt
% 
%t: Zeitpunkte t_k in [t0, T]
%y: Approximationen zu Zeitpunkten t
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 t=t0;
 y=y0;
 N=(T-t0)/h; 
 for i=1:N
  f_=f(t(i),y(:,i));
  y_temp=y(:,i)+h*f_;
  y(:,i+1)=y(:,i)+h/2*(f_+f(t(i)+h,y_temp));
  t(i+1)=t(i)+h;
 end
 end
 
 function[t,y]=richardson_heun(f,t0,y0,h,T,tol)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Verfahren von Heun mit Richardson Schrittweitensteuerung
% zum Loesen der ODE y'=f(t,y), y(t0)=y0
% f: Funktionshandle auf f(t,y)
% t0: Anfangszeit
% y0: Anfangswerte
% h: Schrittweite
% T: Endzeitpunkt
% tol: Fehlertoleranz
%
%t: Zeitpunkte t_k in [t0, T]
%y: Approximationen zu Zeitpunkten t
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 t=t0;
 y=y0;
 ct=1;
 p=%%% % Ordnung des Verfahrens
 Fak=0.8;
 while t(ct)<T
   %Ein Heunschritt mit Schrittweite h
   f_=%%%
   y_temp=%%%
   y(:,ct+1)=%%%
   
   %Zwei Heunschritte mit Schrittweite h/2
   y_temp=%%%
   y_1= %%%
   
   f_=%%%
   y_temp=%%%
   y_tilde= %%%
   
   %Schaetzen des Fehlers
   norm_l=%%%
   
   %Schritt akzeptieren?
   if norm_l<=tol
    t(ct+1)=t(ct)+h;
    ct=ct+1;   
   end
   
   %Vorschlag der Schrittweite
   h=%%%
 end
 end

function[t,y]=trapezmethode(f,t0,y0,h,T,Df)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Trapezmethode zum Loesen der ODE y'=f(t,y), y(t0)=y0
% f: Funktionshandle auf f(t,y)
% y0: Anfangswerte
% t0: Anfangszeit
% h: Schrittweite
% T: Endzeitpunkt
% Df: Funktionshandle auf Jacobimatrix von f
%t: Zeitpunkte t_k in [t0, T]
%y: Approximationen zu Zeitpunkten t
%
% Bem: Das Gleichungsystem der (impliziten) Trapezmethode
% wird mit dem Newtonverfahren geloest mit einem expliziten
% Eulerschritt als Startwert.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 t=t0;
 y=y0;
 tol=1e-6;
 maxit=100;
 N=(T-t0)/h;
 for i=1:N 
  %Newtonfunktion
  fct=@(x) y(:,i)+h/2*(f(t,y(:,i))+f(t+h,x))-x;
  %Jacobimatrix
  Dfct=@(t,y) -eye(length(y0))+h/2*Df(t,y);
  
  %Anfangsbed. fuer Newton (Expl. Euler)
  y_alt=y(:,i)+h*f(t(i),y(:,i));
  
  for k=1:maxit
   %Newtonschritt
   delta=-Dfct(t(i)+h,y_alt)\fct(y_alt);
   y_neu=y_alt+delta;
   if norm(y_alt-y_neu)<tol*(1+norm(y_neu))
    y(:,i+1)=y_neu;
    break;
   end
   y_alt=y_neu;
  end
  t(i+1)=t(i)+h;
 end
end
