%% Numerik D-MAVT, RW (Nipp), SS06
%
% This is the M-file solution to the tenth exercise sheet.
% To run, save this file to your account and ensure that the "Current directory"
% in MATLAB is the same as that where the file is saved.
%
%
% Then to run the file type "s9" in the command window


function s9
% Aufgabe 1
disp(['Aufgabe 1'])
format long e;
x=pi/4;
exact=dfexact(x);
fp_er=[];
for k=1:14
  h=(1/10)^k;  
  fp=(f(x+h)-f(x-h))/(2*h);
  er=abs(exact-fp)/exact;
  fp_er=[fp_er; fp er];
end
fp_er

%----------------------------------------------------------------
% output
% exact =
%      4.444444444444444e-01
% fp_er =
%      4.434320732718755e-01     2.277835138280060e-03
%      4.444345654212088e-01     2.222780228003374e-05
%      4.444443456788071e-01     2.222226840042918e-07
%      4.444444434564818e-01     2.222915851390894e-09
%      4.444444444307025e-01     3.091939204669103e-11
%      4.444444444695603e-01     5.651067114254005e-11
%      4.444444440254710e-01     9.426901759201911e-10
%      4.444444534623667e-01     2.029032504513584e-08
%      4.444444812179423e-01     8.274037018030089e-08
%      4.444444812179422e-01     8.274036993050071e-08
%      4.444444812179422e-01     8.274036993050071e-08
%      4.444222767574496e-01     4.987729573845134e-05
%      4.440892098500620e-01     7.992778373604320e-04
%      4.440892098500620e-01     7.992778373604320e-04
%------------------------------------------------------------------

disp(' Aufgabe 2')
% Aufgabe 2 a
for tf=1:1:3
  x = exp(-0.5*tf); % exact
  disp(['tf=',num2str(tf),', exact=',num2str(x)]);
  for h=[0.1 0.01 0.001] 
    xeuler=eulerexp('fa',tf,h,1);
    disp([' h=',num2str(h),': xeuler=',num2str(xeuler),...
	  ', rel err=',num2str(abs((xeuler-x)/x))]);
  end
end
x=-0.23256716;	% b)exact
tf=1;
disp(['b): exact=',num2str(x)]);
for h=[0.1 0.01 0.001] 
  xeuler=eulerexp('fb',tf,h,-1);
  disp([' h=',num2str(h),': xeuler=',num2str(xeuler),...
        ', rel err=',num2str(abs((xeuler-x)/x))]);
end

%-----------------------------------------------------------------
% output
% a):
% tf=1, exact=0.60653
%  h=0.1: xeuler=0.59874, rel err=0.01285
%  h=0.01: xeuler=0.60577, rel err=0.0012534
%  h=0.001: xeuler=0.60645, rel err=0.00012503
% tf=2, exact=0.36788
%  h=0.1: xeuler=0.35849, rel err=0.025534
%  h=0.01: xeuler=0.36696, rel err=0.0025052
%  h=0.001: xeuler=0.36779, rel err=0.00025005
% tf=3, exact=0.22313
%  h=0.1: xeuler=0.21464, rel err=0.038056
%  h=0.01: xeuler=0.22229, rel err=0.0037555
%  h=0.001: xeuler=0.22305, rel err=0.00037505
% b): exact=-0.23257
%  h=0.1: xeuler=-0.24285, rel err=0.044225
%  h=0.01: xeuler=-0.2337, rel err=0.0048769
%  h=0.001: xeuler=-0.23268, rel err=0.00049211
%----------------------------------------------------------------

disp(['Aufgabe 4 working...this may take some time:-)'])
y0=[0.5;0;0;sqrt(3)];tf=500; 	% initial condition, final time
% compute using ode23
[T,Y]=ode23('twobody',[0 tf],y0);
% plot resulting orbit
plot(Y(:,1),Y(:,2)); title('using ode23');
% plot the stepsizes
figure;
plot(diff(T),'.-'); title('stepsizes h from ode23');
xlabel('i'); ylabel('h_i');

% compute using Heun
for h=[0.1 0.01]
  % integrate using heun
  y = y0; orbit=y(1:2);
  for t=0:h:tf-h
    y = heun(t,y,h,'twobody');
    orbit = [orbit y(1:2)];
  end
  % plot orbit and trajectories
  figure;  
  plot(orbit(1,:),orbit(2,:),'.:');
  xlabel('y_1(t)');ylabel('y_2(t)');
  title(['h=',num2str(h)]);
  figure;  
  plot([0:h:tf],orbit(1,:),':',[0:h:tf],orbit(2,:),'-.');
  legend('y_1(t)','y_2(t)');
  xlabel('t');
  title(['h=',num2str(h)]);
end

%------------------------------------------------------------------
function ff=f(x,a,b)
% function for 1)
if nargin==1 a=1;b=0; end
ff=1./(1+cos(a*x+b).^2);
%------------------------------------------------------------------
function df=dfexact(x,a,b)
% exact derivative function for 1)
if nargin==1 a=1;b=0; end
df=2*a*cos(a*x+b).*sin(a*x+b).*f(x,a,b).^2;
%------------------------------------------------------------------
function xf = eulerexp(f,tf,h,x0)
% explicit euler for d/dt x(t) = f(t,x), x(0)=x0, step h
x = x0; 
for t=0:h:tf-h,
  x = x + h*feval(f,t,x);
end
xf = x;
%------------------------------------------------------------------
function f=fa(t,x)
% function for 2a)
f=-0.5*x;
%------------------------------------------------------------------
function f=fb(t,x)
% function for 2b)
f=t.^2+x.^2;
%-------------------------------------------------------------------
function ykp1 = heun(tk,yk,h,f)
% one step of Heun procedure with stepsize h:
% ykp1 = yk + h/2(f(tk,yk)+f(tk+h,yk+h*f(tk,yk)))
fk   = feval(f,tk,yk);
ykp1 = yk + h/2*(fk+feval(f,tk+h,yk+h*fk));
%-------------------------------------------------------------------
function f=twobody(t,y)
% function for 4
d=(y(1)*y(1)+y(2)*y(2))^(1.5);
f=[y(3);y(4);-y(1)/d;-y(2)/d];
