Midpoint Method on an SHO to a Velocity Verlet Method

t(1,1) = 0;    % Starting time (arbitrary units)
u(1,1) = 1;    % Initial position (arbitrary units)
u(1,2) = 0;    % Initial velocity (arbitrary units)
tend = 4*pi; % Simulation run time.
Nsteps = 200;  % Number of integration steps;
dt = (tend-t(1))/Nsteps;  % Time step size (arbitrary units)

%% Set up the differential equation
% This is the "right hand side" function of the differential equation, T'
k = 1; m = 1; omega = sqrt(k/m);       % Physical parameters of oscillator
dudt = @(u) [ u(2) -omega^2 * u(1) ];  % Right-hand side function

% Midpoint Method

%for idx = 1:tend/dt                      
%    t(idx+1,1) = t(idx,1) + dt;        
%    utmp = u(idx,:) + (dt/2)*dudt(u(idx,:));
%    u(idx+1,:) = u(idx,:) + dt*dudt(utmp);
%end

%Velocity Verlet Method

for idx = 1:tend/dt
    t(idx+1,1) = t(idx,1) + dt;
    uv = u(:,idx) - (dt/2)*dudt(u(idx,:));
    u(idx+1,:) = u(idx,:) + dt*dudt(uv);
    u(:,idx+1) = uv - (dt/2)*dudt(u(idx+1,:));
end

Here is my code I am try to fix. I have the midpoint method for a simple harmonic oscillator, which I know runs and works. I am curently trying to edit it and make it into a velocity verlet method but I can not figure out what I have done wrong, and if I am not even close please let me know as well. Any and all help will be appreciated.

Answer 1

[Dongyue]

The following commands will lead to dimension inconsistency. Could you please verify whether they are correct? You may exchange row indexing and column indexing for the matrix u.

u(:,idx) - (dt/2)*dudt(u(idx,:));

u(:,idx+1) = uv - (dt/2)*dudt(u(idx+1,:));