clear all; close all;


    % Set configuration parameters
    mass = 1;           % Particle mass
    dt = 0.0001;        % Integration time
    dt2 = dt*dt;        % Integration time, squared
    dwPar = 1;          % A parameter to scale the double well potential

    % Set simulation parameters
    nSteps = 100000;    % Total simulation time (in integration steps)
    sampleFreq = 1;     % Sampling frequency
    sampleCounter = 1;  % Sampling counter

    % Set trajectories to follow particle's position and velocity
    xTraj = zeros(1,nSteps/sampleFreq);
    vTraj = zeros(1,nSteps/sampleFreq);

    % Set initial conditions
    x = 0;              % Initial position
    oldX = -0.001;      % Position in previous time step


    % ===================
    % Molecular Dynamics
    % ===================
    for step = 1:nSteps

        % Tmp will hold x(t)
        tmp = x;

        % Calculate the new x positon by integrating the equations of motion
        % x(t+dt) = 2*x(t) - x(t-dt) + dt^2*(f(t)/m) + O(dt^4)
        x = 2.0*x - oldX + (-1)*dwPar*4*x.*(x.^2-1)*dt2/mass;

        % oldX now holds x(t), and x holds x(t+dt)
        oldX = tmp;

        % Integrate the velocity as well.
        % O(dt) accurate. See Verlet integration on wikipedia
        % v(t+dt) = (x(t+dt) - x(t))/ dt + O(dt)
        v = (x - oldX)/dt;

        % Sample
        if mod(step,sampleFreq) == 0
            xTraj(sampleCounter) = x;
            vTraj(sampleCounter) = v;
            sampleCounter = sampleCounter + 1;
        end

    end

    % ===================
    % Simulation results
    % ===================

    % Plot the particle's position trajectory
    figure(1)
    plot(xTraj);
    xlabel('time')
    ylabel('Position')

    % Plot the sum of kinteic and potential energy - to check for conservation
    figure(2)
    potEnergy = dwPar*(xTraj.^2-1).^2;
    kinEnergy = 0.5*mass*(vTraj.^2);
    plot(potEnergy + kinEnergy);
    xlabel('time')
    ylabel('Total Energy')

Published with MATLAB® R2012b