Problem with solving discrete element method using leap frog method

Thin Rupar Win · Jan 15, 2024 · 2.3K views
Question
I am new to writing matlab programming by using Discrete element method using leap frog algorithm. I got many error coming from my program. Can you all suggest me how to correct them with your all idea? Please let me hear your reply.    n_part=4; kn=5; kt=2/7*kn; m=0.3; g=9.81; rad(1:n_part)=0.5; v_init=0.5; y=zeros(); % testing i_particle=1:n_part; % testing j_particle=n_part+1:2*n_part; position_x(1,1:n_part)=0.05; theta=2*pi*rand(size(position_x)); velocity_x=v_init*sin(theta); position_y(1,1:n_part)=0.05; velocity_y=v_init*cos(theta); timestep=100; dt=0.001; acceleration_x(:,1:n_part)=zeros(); acceleration_y(:,1:n_part)=zeros(); Fn=zeros(); Fn_i=zeros(); Fn_j=zeros(); v_half_x(1,1:n_part)=zeros(); v_half_y(1,1:n_part)=zeros(); for n=1:n_part for k=2:timestep % position_x v_half_x(k,n)=velocity_x(k-1,n)+0.5*dt*acceleration_x(k-1,n); position_x(k,n)=position_x(k-1,n)+v_half_x(k-1,n)*dt; % position_y v_half_y(k,n)=velocity_y(k-1,n)+0.5*dt*acceleration_y(k-1,n); position_y(k,n)=position_y(k-1,n)+v_half_y(k-1,n)*dt; for i=1:n_part for j=i+1:n_part if i>j % real position & distance lx=position_x(k-1,i)-position_x(k-1,j); ly=position_y(k-1,i)-position_y(k-1,j); root_xy=sqrt(ly^2+ly^2); % force calculation Fn=kn*root_xy^1.5; Fn_i=Fn_i+Fn; Fn_j=Fn_j+Fn; % acceleration term acceleration_x(k,:)=Fn_i./m; acceleration_y(k,:)=Fn_j./m; end end end velocity_y(k,n)=v_half_y(k-1,n)+0.5*dt*acceleration_y(k-1,n); velocity_x(k,n)=v_half_x(k-1,n)+0.5*dt*acceleration_x(k-1,n); end end  
Expert Answer
Profile picture of John Williams
John Williams PhD Expert
Answered Aug 23, 2026

The following gets the code working, but I've no idea if the results are meaningful!!

 

n_part=4;
kn=5;
kt=2/7*kn;
m=0.3;
g=9.81;
rad(1:n_part)=0.5;
v_init=0.5;

y=0;
% testing i_particle=1:n_part;
% testing j_particle=n_part+1:2*n_part;

position_x(1,1:n_part)=0.05;
theta=2*pi*rand(size(position_x));
velocity_x=v_init*sin(theta);
position_y(1,1:n_part)=0.05;
velocity_y=v_init*cos(theta);
timestep=100;
dt=0.001;
acceleration_x=zeros(timestep,n_part); %%%%%%%%%%%%%%
acceleration_y=zeros(timestep,n_part); %%%%%%%%%%%%%%
Fn=0;
Fn_i=0;
Fn_j=0;
v_half_x=zeros(timestep,n_part); %%%%%%%%%%%%%%
v_half_y=zeros(timestep,n_part); %%%%%%%%%%%%%%
for n=1:n_part
    for k=2:timestep
        % position_x
        v_half_x(k,n)=velocity_x(k-1,n)+0.5*dt*acceleration_x(k-1,n);
        position_x(k,n)=position_x(k-1,n)+v_half_x(k-1,n)*dt;
        % position_y
        v_half_y(k,n)=velocity_y(k-1,n)+0.5*dt*acceleration_y(k-1,n);
        position_y(k,n)=position_y(k-1,n)+v_half_y(k-1,n)*dt;
        
        for i=1:n_part
            for j=i+1:n_part
                %if i>j    %%%%% i CANNOT be greater than j as you %%%%%%
                           %%%%% set j to be i+1 upwards!          %%%%%%
                    lx=position_x(k-1,i)-position_x(k-1,j);
                    ly=position_y(k-1,i)-position_y(k-1,j);
                    root_xy=sqrt(ly^2+ly^2);
                    % force calculation
                    Fn=kn*root_xy^1.5;
                    Fn_i=Fn_i+Fn;
                    Fn_j=Fn_j+Fn;
                
                    % acceleration term
                    acceleration_x(k,:)=Fn_i./m;
                    acceleration_y(k,:)=Fn_j./m;
               % end
            end
        end
        velocity_y(k,n)=v_half_y(k-1,n)+0.5*dt*acceleration_y(k-1,n);
        velocity_x(k,n)=v_half_x(k-1,n)+0.5*dt*acceleration_x(k-1,n);
    end
end

 

100% Run Guarantee 3-Hour Fast-Track Delivery

Need a Custom Version or Complete Simulation for This Problem?

Our 500+ PhD engineers build, debug, and optimize working MATLAB scripts and Simulink (.slx) models tailored to your exact assignment rubrics with zero plagiarism.

Tested on MATLAB R2024b / R2026a
Turnitin 0% Plagiarism Report
Free 7-Day Revisions Guarantee
Have a different question? Ask here

Get a Free Consultation or a Sample Assignment Review!