% script: FFTest.m

% Octave and Matlab compatible file
disp('FFT experiments:')
disp('dy/dt = fun(t,y),  0<=t<=tLast')
disp('y(0) = y0')
% initialization code
t0 = 0;
tLast = 1;
N  = 16
D = 0.5;
xNodes = linspace(0,2*pi,N+1)'; % make sure it's a column
xNodes = xNodes(1:N);
global gy0;
gy0 =  sin(xNodes).^2;
fun = @(t,y) -exp(-t/2)/2*(2*D*cos(2*xNodes) + sin(xNodes).^2);
funCap = @(t,y) FFTransform(fun(t,y),0);
numSteps = 30*[1,2,4,8];
hSteps = (tLast - t0)./numSteps;
lenSteps = length(hSteps);
errors = zeros(lenSteps,3);
% store stepsizes and errors in this matrix
errors(:,1) = hSteps;
% main  -- we use the infinity norm for error
ySoln = FFTestSoln(tLast);
for jj = 1:lenSteps
  yApprox = Euler(fun, t0, gy0, hSteps(jj), numSteps(jj));
  errors(jj,2) = norm(yApprox - ySoln,inf);
  yApprox = ERK4(fun, t0, gy0, hSteps(jj), numSteps(jj));
  errors(jj,3) = norm(yApprox - ySoln,inf);
end
disp('Table with columns [time stepsize, Euler error, ERK4 error]:');
format long
disp(errors)
disp('Table with columns of rates of decrease of error [stepsize, Euler, ERK4]');
disp(errors(1:lenSteps-1, 1:3)./errors(2:lenSteps, 1:3))
format

% now repeat the experiment using FFTs
y0Cap = FFTransform(gy0, 0);  % original y0 transformed
% main  -- we use the infinity norm for error
for jj = 1:lenSteps
  yApproxCap = Euler(funCap,t0,y0Cap,hSteps(jj),numSteps(jj));
  yApprox = FFTransform(yApproxCap,1);
  errors(jj,2) = norm(yApprox - ySoln,inf);
  yApproxCap = ERK4(funCap,t0,y0Cap,hSteps(jj),numSteps(jj));
  yApprox = FFTransform(yApproxCap,1);
  errors(jj,3) = norm(yApprox - ySoln,inf);
end
disp(' ');
disp('Repeat experiment with FFTs: ')
disp('Table with columns [time stepsize, Euler error, ERK4 error]:');
format long
disp(errors)
disp('Table with columns of rates of decrease of error [stepsize, Euler, ERK4]');
disp(errors(1:lenSteps-1, 1:3)./errors(2:lenSteps, 1:3))
format


