Midpoint Method on an SHO to a Velocity Verlet Method

6 次查看（过去 30 天）
Drake
Drake 2021-2-21
回答： Dongyue 2022-11-17
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 Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by