% M-File grxxm073.m
% Adams-Bashforth-Verfahren
% y' = - y + cos(at),  y(0) = y_0,    a  reeller Parameter

% Grafikfenster (nur MATLAB / OCTAVE: man kommentiere die naechsten 3 Zeilen)
  screen=get(0,'ScreenSize');
  screenplot=[50 50 screen(3)*0.8 screen(4)*0.8];
  figure('Position',screenplot,'color',[1 1 0.77]);

clf;
clear;

% Koeffizientenmatrix für Adams-Bashforth (Spalten vertauscht)
K=[0,0,0,-1/2,3/2
     0,0,5/12,-16/12,23/12
     0,-9/24,37/24,-59/24,55/24
     251/720,-1274/720,2616/720,-2774/720,1901/720];

% Differentialgleichung
f=inline('-y+cos(a*t)','t','y','a');

% Loesung
l=inline('(y0-1/(1+a^2))*exp(-t)+(cos(a*t)+a*sin(a*t))/(1+a^2)','t','y0','a');

% Aufgabendaten
y0=2;
a=8;
T=2*pi;

% Schrittweite h
n=200;     %%   Verfahrensparameter  %%
h=T/n; 

% Ausgabevorbereitung
hnum=num2str(h);
anum=num2str(a);

% r=1 
% Anlaufrechnung mit Heun-Verfahren
y1(1,1:5)=[0 y0 -y0+1 y0 0];
y1(2,1)=h;
    k1=f(y1(1,1),y1(1,2),a);
    k2=f(y1(1,1)+h,y1(1,2)+h*k1,a);
y1(2,2)=y1(1,2)+h*(k1+k2)/2;
y1(2,3)=f(y1(2,1),y1(2,2),a);
y1(2,4)=l(y1(2,1),y0,a);
y1(2,5)=abs(y1(2,2)-y1(2,4));
for k=3:n
    y1(k,1)=k*h;
    y1(k,2)=y1(k-1,2)+h*K(1,4:5)*y1((k-2):(k-1),3);
    y1(k,3)=f(y1(k,1),y1(k,2),a);
    y1(k,4)=l(y1(k,1),y0,a);
    y1(k,5)=abs(y1(k,2)-y1(k,4));
end;
subplot(2,1,1);
plot(y1(:,1),y1(:,4),'b',y1(:,1),y1(:,2),'r',[0 T],[0 0],'k');
legend('exakte Lösung','A-B-Verfahren');
str1='Adams-Bashforth-Verfahren : ';
str2=[' r = 1,  a = ' anum  ' , h = ' hnum  ',    Anlaufrechnung:  Heun-Verfahren'];
title([str1 str2],'FontSize',15);


% r=2 
% Anlaufrechnung mit Runge-Kutta-Verfahren
y2(1,1:5)=[0 y0 -y0+1 y0 0];
for s=1:2
    y2(s+1,1)=s*h;
       k1=f(y2(s,1),y2(s,2),a);
       k2=f(y2(s,1)+h/2,y2(s,2)+h*k1/2,a);
       k3=f(y2(s,1)+h/2,y2(s,2)+h*k2/2,a);
       k4=f(y2(s,1)+h,y2(s,2)+h*k3,a);
    y2(s+1,2)=y2(s,2)+h*(k1+2*k2+2*k3+k4)/6;  
    y2(s+1,3)=f(y2(s+1,1),y2(s+1,2),a);
    y2(s+1,4)=l(y2(s+1,1),y0,a);
    y2(s+1,5)=abs(y2(s+1,2)-y2(s+1,4));
end;
for k=4:n
    y2(k,1)=k*h;
    y2(k,2)=y2(k-1,2)+h*K(2,3:5)*y2((k-3):(k-1),3);
    y2(k,3)=f(y2(k,1),y2(k,2),a);
    y2(k,4)=l(y2(k,1),y0,a);
    y2(k,5)=abs(y2(k,2)-y2(k,4));
end;
subplot(2,1,2);
plot(y2(:,1),y2(:,4),'b',y2(:,1),y2(:,2),'r',[0 T],[0 0],'k');
legend('exakte Lösung','A-B-Verfahren');
str3=[' r = 2,  a = ' anum  ' , h = ' hnum  ',    Anlaufrechnung:  R-K-verfahren'];
title([str1 str3],'FontSize',15);

% Ausgabe des absoluten Fehlers
% >> format long
% >> [y1(:,5) y2(:,5)]
