Problem with solving discrete element method using leap frog method

Dear Sir or Madam,

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

Answer 1

[Alan Stevens]

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

[Thin Rupar Win]

Thank you very much for your idea. I acknowledge your help. Have a nice day.