% testspring  -  test numerical solvers using simple mass-spring model

% The differential equation is   x" + x = 0.
% The time interval of interest is [0,tMax], with tMax chosen near the top
% The initial values x(0)=x0 and x'(0)=v0 are defined near the top.

% (c) 2014, Philip D Loewen; history at end of file

% Define time interval of interest and initial state
tMax   = 8*pi;
x0     = 1;
v0     = 1;

% Defining y=[x;v] turns x"+x=0 into the system y'=G(t,y) for this G:
G = @(t,y)  [y(2); -y(1)];
y0 = [x0;v0];     % The initial information fills this vector.

star16 = '****************';
ode45options = odeset('RelTol',1.e-7,'AbsTol',1.e-9);

% LOOP: Run through each of the selected methods
%for mmm = [1, 2, 3]   % mmm sets what method to use. Details below.
for mmm = [3]
  switch mmm
    case 1
      method = 'Euler';
    case 2
      method = 'Heun';
    case 3
        method = 'Runge-Kutta plus (ode45)';
    otherwise
        method = 'Runge-Kutta plus (ode45)';  % our fallback default
  end
  
  disp([star16,' New Method: ',method,' ',star16]);
 
  % LOOP:  Apply the selected method to each mesh size specified next
  for Nsubint=[10, 20]
    tspan = linspace(0,tMax,1+Nsubint);

    switch mmm
      case 1
        [Tcol,Ycol] = ode01(G,tspan,y0);   % Students must supply ode01 [Euler]
      case 2
        [Tcol,Ycol] = ode12(G,tspan,y0);   % Students must supply ode12 [Heun]
      case 3
        [Tcol,Ycol] = ode45(G,tspan,y0,ode45options);
      otherwise
        [Tcol,Ycol] = ode45(G,tspan,y0);
    end

    % Any students not named Potter should put his/her own name into the title.
    titletext = sprintf([method,' Method with %d Subintervals [Potter]'],...
                    Nsubint);

    figure('Name',titletext);
    annotation('textbox', [0 0.9 1 0.1], ...
      'String', titletext, ...
      'EdgeColor', 'none', ...
      'HorizontalAlignment', 'center');

    vm = 0.1; % vertical margin
    hm = 0.1; % horizontal margin

    phaseplane = subplot('position',[0.0+hm,0.0+vm,0.6-2*hm,1.0-2*vm]); hold on;
    timetraj   = subplot('position',[0.6+hm,0.5+vm,0.4-2*hm,0.5-2*vm]); hold on;
    energy     = subplot('position',[0.6+hm,0.0+vm,0.4-2*hm,0.5-2*vm]); hold on;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    axes(phaseplane);
    axis equal;
    set(gca,'XLim',2*[-1,1],'YLim',2*[-1,1]);
    grid on;
    hold on;
    title('Phase Plane View of Mass-Spring State');
    xlabel('Spring Displacement x, m');
    ylabel('Spring Velocity v, m/s');

    plot(Ycol(:,1),Ycol(:,2),'r','LineWidth',1);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    axes(timetraj);
    axis equal;
    set(gca,'XLim',tMax*[0,1],'YLim',2*[-1,1]);
    grid on;
    hold on;
    title('Mass position and velocity vs time');
    xlabel('Time t, s');
    ylabel('Position x, m, and Velocity v, m/s');

    plot(Tcol,Ycol(:,1),'r');
    plot(Tcol,Ycol(:,2),'b');
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    axes(energy);
   
    grid on;
    xlabel('Time t, s');
    ylabel('Energy E = (x^2+v^2)/2, J');
    title('Total Energy [T+V]');
    set(gca,'XLim',tMax*[0,1]);

    E = (Ycol(:,1).^2 + Ycol(:,2).^2)/2;
    plot(Tcol,E);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    disp(sprintf('Subinterval count: %4d',Nsubint));
    disp(sprintf('Final time t:      %13.8f',tMax));
    disp(sprintf('Final position x:  %13.8f (true sol: %13.8f)',...
                 [Ycol(1+Nsubint,1),x0*cos(tMax)+v0*sin(tMax)]));
    disp(sprintf('Final velocity v:  %13.8f (true sol: %13.8f)',...
                 [Ycol(1+Nsubint,2),-x0*sin(tMax)+v0*cos(tMax)]));
    disp(sprintf('Final energy T+V:  %13.8f (true sol: %13.8f)',...
                 [E(1+Nsubint),(x0^2+v0^2)/2]));
    disp(' ');
  end
end;
    
% Version history (reverse order)
%  2014-11-22: First draft, adapted from newpendulum.m
