# MATLAB Piecewise MATLAB ode45 problem for Lennard-Jones

1. Feb 16, 2015

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')

#### Attached Files:

File size:
37.6 KB
Views:
325
• ###### Chapter 8 Hard Boundary.docx
File size:
22.7 KB
Views:
110
2. Feb 21, 2015

### Greg Bernhardt

Thanks for the post! This is an automated courtesy bump. Sorry you aren't generating responses at the moment. Do you have any further information, come to any new conclusions or is it possible to reword the post?

3. Feb 23, 2015

Do you mean that the post cant be comprehended?

4. Feb 23, 2015

### Staff: Mentor

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).

5. Feb 27, 2015

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

6. Feb 27, 2015

### Staff: Mentor

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.