Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Piecewise MATLAB ode45 problem for Lennard-Jones

  1. Feb 16, 2015 #1
    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:

  2. jcsd
  3. Feb 21, 2015 #2
    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?
     
  4. Feb 23, 2015 #3
    Do you mean that the post cant be comprehended?
     
  5. Feb 23, 2015 #4

    DrClaude

    User Avatar

    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).
     
  6. Feb 27, 2015 #5
    Basically, t = 4.1 s means that the particle collides with the hard boundary by using ode45.
     
  7. Feb 27, 2015 #6

    DrClaude

    User Avatar

    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Piecewise MATLAB ode45 problem for Lennard-Jones
  1. Matlab ODE45 Help? (Replies: 4)

  2. Matlab ode45 (Replies: 7)

  3. MATLAB help-ode45 (Replies: 0)

Loading...