Piecewise MATLAB ode45 problem for Lennard-Jones

[Adrian Soo]

I could run the MATLAB codes but the main issue is that data for T1 and U1 could not be produced. I am currently trying to calculate the displacement of two particles after the collision with walls with length of 10 units.

Here are the codes. If you don't understand its context, you can refer to the Microsoft word files I posted. Note that each particle HAS A RADIUS OF 1 UNIT.​

================================================================================
Lennard-Jones's equation of motion function

function du = lennardJones(~,u)

du = zeros(8,1);

du(1) = u(2);
du(2) = -24*( u(1) - u(3) )/( ( u(1) - u(3) )^2 + ( u(5) - u(7) )^2 )^7 + 24*( u(1) - u(3) )/( ( u(1) - u(3) )^2 + ( u(5) - u(7) )^2 )^4;
du(3) = u(4);
du(4) = -24*( -u(1) + u(3) )/( ( u(1) - u(3) )^2 + ( u(5) - u(7) )^2 )^7 + 24*( -u(1) + u(3) )/( ( u(1) - u(3) )^2 + ( u(5) - u(7) )^2 )^4;
du(5) = u(6);
du(6) = -24*( u(5) - u(7) )/( ( u(1) - u(3) )^2 + ( u(5) - u(7) )^2 )^7 + 24*( u(5) - u(7) )/( ( u(1) - u(3) )^2 + ( u(5) - u(7) )^2 )^4;
du(7) = u(8);
du(8) = -24*( -u(5) + u(7) )/( ( u(1) - u(3) )^2 + ( -u(5) + u(7) )^2 )^7 + 24*( -u(5) + u(7) )/( ( u(1) - u(3) )^2 + ( u(5) - u(7) )^2 )^4;================================================================================
clear all

t0 = 0; % Initial time
deltat = 0.1; % Step size of time
tend = 10; % Final time
u = [4,1,6,-1,1,0,1,0]; % Initial conditions for Lennard-Jones Potential
n = (tend - t0)/deltat; % for 'for' loops
t1 = t0; % Initialise t1 to t0
r = 1;

for i = 1 : n - 1

t2 = t1 + deltat; % Increment in time
options = odeset('RelTol',1e-8,'AbsTol',[1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8 1e-8]);
[T,U] = ode45(@lennardJones,[t1 t2],u,options); % ode45 applied on LennardJones


if U(1) <= 1 && U(3) >= 9 %When particle collides with the wall at the distance of 9 units on the right
% and 1 unit on the left
clear t1 t2 deltat
t1 = 4.1; % t1 indicates the time it collides with the solid wall
deltat = 0.1; % increment in time
t2 = 10;
k = (t2 - t1)/deltat; %for 'for' loops
u1 = [1,1,9,-1,1,0,1,0]; %initial condition changes after colliding with walls
newtime = 0;

for j = 1:k
newtime = newtime + deltat;
[T1,U1] = ode45(@lennardJones,[t1 newtime],u1,options);
end

end

end

plot(T,U,T1,U1,'o')
============================================================================

output :
Undefined function or variable 'T1'.

Error in hardBoundary (line 37)
plot(T,U,T1,U1,'o')

 

[Adrian Soo]

Do you mean that the post can't be comprehended?

 

[DrClaude]

This is an automated courtesy bump.

This message is sent automatically to new posts that do not have any replies after one week.

I don't understand how your program is supposed to deal with collisions. It seems to be integrating over a certain time interval, and then come back to an arbitrary point in time (t=4.1) if there is a collision. Could you describe in more details?

Also, please use code tags when posting code: [C O D E] Some code here [/C O D E] (without spaces in CODE).

 

[Adrian Soo]

Basically, t = 4.1 s means that the particle collides with the hard boundary by using ode45.

 

[DrClaude]

Basically, t = 4.1 s means that the particle collides with the hard boundary by using ode45.

How do you know?

As for your original problem: the condition
if U(1) <= 1 && U(3) >= 9
is never met, hence T1 is never defined.