Strange ODE solution to Damped Driven Harmonic Oscillator

Problem Description:

Im looking into force damped harmonic oscillation with forcing taking the form of a square wave. I have implemented one basic ode solver myself (see section 3.2 in code) and I am testing this against the solver ode45. What I cannot seem to understand is the phase of the oscillation with respect to the forcing function. It is different in both cases. I am unsure if this is to do with my implemented solver or my implementation with the ode45 solver.

The Code The code is as follows it should be fully functional and includes the two figures I have been using to view the results.

if true
FE=4000;                    % Frequency of Forcing (Hz)
F0=4000;                    % Oscillation Frequency
sigma=700;                  % Damping Factor
Fs=100000;                  % Sampling Frequency
An_Amp=1;                   % Amplitude of Forcing
SirenTeeth=8;               % Number of teeth before decay
MissingTeeth=7;             % Number of teeth not there during decay
Inital_Amp=0;               % Inital amplitude of oscillation
%--------------------------------------------------------------------------
%   1.1 Analysis Parameters
%--------------------------------------------------------------------------
tstart=0;                   % Time Start (sec)
tend=0.1;                   % Time End (sec)
% Resolution needs to be higher than the sampling frequency we use. 
Calc_res=10;               % Increases the resolution of the calculation by the factor
%--------------------------------------------------------------------------
%   1.1 Calculated Parameters
%--------------------------------------------------------------------------
w0=F0*2*pi;                 % Angular frequency of system
gamma=2*sigma;              % Full Width at Half Maximum
dt=1/(Fs*Calc_res);         % Time step
Duration=tend-tstart;       % Duration of investigation
Time=(tstart:dt:tend)';     % Time vector
SirenPeriod=(SirenTeeth+MissingTeeth)/FE;
%__________________________________________________________________________

2. Data Pre Processing

%--------------------------------------------------------------------------
%   2.1 Prepare for interations
%--------------------------------------------------------------------------
x1(1)=Inital_Amp;            % Starting oscillation amplitude
xp=0;                        % Starting gradent
%__________________________________________________________________________

3. Data Prosessing

%--------------------------------------------------------------------------
%   3.1 Driving Function
%--------------------------------------------------------------------------
%possible issues with this method for quickly ramping excitation
%frequencies
for n=1:(SirenPeriod/dt);
  if n<(SirenPeriod/dt)*(SirenTeeth/(SirenTeeth+MissingTeeth))
  Anregung(n)=An_Amp*round(0.5+0.5*sin(2*pi*FE*Time(n:n)));
  antiregung(n)=0;
  else
  Anregung(n)=0;
  antiregung(n)=An_Amp*round(0.5+0.5*sin(2*pi*FE*Time(n:n)));
  end
end
Anregung2=[];
antiregung2=[];
for n=1:(SirenPeriod/dt):(length(Time))
  Anregung2=[Anregung2,Anregung];
  antiregung2=[antiregung2,antiregung];
end
Anregung=Anregung2(:,1:length(Time));
Antiregung=antiregung2(:,1:length(Time));
%--------------------------------------------------------------------------
%   3.2(forward trapp?)
%--------------------------------------------------------------------------
for n=1:(length(Time)-1)
  xpp=-gamma*xp-w0^2*x1(n)+Anregung(n);
  dxp=xpp*dt;
  dx=xp*dt;
  xp=xp+dxp;
  x1(n+1)=x1(n)+dx;
end
x1=x1';
%--------------------------------------------------------------------------
%   3.4 Rungekutta matlab
%--------------------------------------------------------------------------    
  %----------------------------------------------------------------------
  %   3.4.1 Initial conditions
  %----------------------------------------------------------------------
  x=0;
  xp=0;
  %----------------------------------------------------------------------
  %   3.4.2 Linearise equations
  %----------------------------------------------------------------------
  function [Ydot]=Oscl(t,Y,gamma,Anregung,w0,Time)    
  Anregung_now=interp1(Time,Anregung,t,'nearest'); 
  Ydot=zeros(2,1);
  Ydot(1)=Y(2);
  Ydot(2)=-gamma*Y(2)-w0^2*Y(1)+Anregung_now;

    end
    %----------------------------------------------------------------------
    %   3.1 Solve ODE
    %----------------------------------------------------------------------
    % (t,Y) are there to specifiy the order of input arguments which is
    % required when the func
    [T,Y]=ode45(@(t,Y) Oscl(t,Y,gamma,Anregung,w0,Time),Time',[x xp]);
%__________________________________________________________________________

4. Data Output/Save Data Output

%--------------------------------------------------------------------------
%   4.1 Plotting
%--------------------------------------------------------------------------
scrsz=get(0,'ScreenSize');
Freq_min=0.8*min(FE);
Freq_max=1.2*max(FE);
Time_cent_start=tstart+0.5*(Duration)-SirenPeriod;
Time_cent_end=tstart+0.5*(Duration)+SirenPeriod;
    %----------------------------------------------------------------------
    %   4.1.1 Amplitude Time with forcing square wave
    %----------------------------------------------------------------------
    F1=figure('Position',[10 (scrsz(4)/2-20) scrsz(3)/2-20 (scrsz(4)/2-60)]);
    F1S1=subplot(1,1,1,'Parent',F1);
    plot(F1S1,Time,x1/max(x1),'DisplayName','Damped Oscillator')
    hold on
    plot(F1S1,Time,Anregung,'DisplayName','Damped Oscillator','Color','r');
    axis([Time_cent_start Time_cent_end -1.1 1.1]);
    xticks = get(F1S1,'XTick');
    set(F1S1,'XTickLabel', sprintf('%.4f|',xticks))
    grid(F1S1,'on')
    legend(F1S1,'Pressure Amplitude','Driving Function','Location','NorthEast')

    Fig1_title_pt1='Amplitude(t) Trap Integration method';
    Fig1_title_pt2=['alpha=',num2str(sigma),'(1/s)','  ',...
                'Driving Frequency=',num2str(min(FE)),'-',num2str(max(FE)),'Hz'];
    title({Fig1_title_pt1; Fig1_title_pt2});
    xlabel('time (s)')
    ylabel('x(t)')
    %----------------------------------------------------------------------
    %   4.1.1 Amplitude Time with forcing square wave
    %----------------------------------------------------------------------
    F2=figure('Position',[scrsz(3)/2+10 (scrsz(4)/2-20) scrsz(3)/2-20 (scrsz(4)/2-60)]);
    F2S1=subplot(1,1,1,'Parent',F2);
    plot(F2S1,T,Y/max(Y),'DisplayName','Damped Oscillator')
    hold on
    plot(F2S1,Time,Anregung,'DisplayName','Damped Oscillator','Color','r')
    axis([Time_cent_start Time_cent_end -1.1 1.1]);
    xticks = get(F2S1,'XTick');
    set(F2S1,'XTickLabel', sprintf('%.4f|',xticks))
    grid(F2S1,'on')
    legend(F2S1,'Pressure Amplitude','Driving Function','Location','NorthEast')

    Fig2_title_pt1='Amplitude(t) RungeKutta method';
    Fig2_title_pt2=['alpha=',num2str(sigma),'(1/s)','  ',...
                'Driving Frequency=',num2str(min(FE)),'-',num2str(max(FE)),'Hz'];
    title({Fig2_title_pt1; Fig2_title_pt2});
    xlabel('time (s)')
    ylabel('x(t)')
  end

Naturally any one who could shed some light on this from a mathematical point of view or highlight the bad parts in my code would make me very happy!

Cheers, Sam

Answer 1

[Sam]

Question Update:

By updating the code to use:

if true
  sol=ode45(@(t,Y) Oscl(t,Y,gamma,Anregung,w0,Time),Time',[x xp]);
  Eval=deval(sol,Time',1);
  plot(Time',Eval/max(Eval))
end

The phase of the matlab ode method now matches that of section 3.2 result. Therefore, im concluding it was an issue with the time vectored solution to the ode45. Any further comments or suggestions are still more than welcome, maybe it will help my understanding.

Cheers, Sam