Midpoint Method on an SHO to a Velocity Verlet Method

3 次查看(过去 30 天)
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.

回答(1 个)

Dongyue
Dongyue 2022-11-17
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,:));

类别

Help CenterFile Exchange 中查找有关 Programming 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by