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

[MATLAB] Modeling response of a particle in an anharmonic potential

  1. Apr 8, 2015 #1
    I'm trying to model the response of a particle in an anharmonic potential in MATLAB. I know that writing position, x as a function of potential energy U is,


    where T is the period, m is mass, and E is the energy of the system. How would I get a function for T(E) since it's not independent of amplitude obviously. As a starting point, I'm using the anharmonic potential of

    Where a and b are constants. Does anybody know how I could program this in matlab?

  2. jcsd
  3. Apr 9, 2015 #2


    User Avatar
    Gold Member

  4. Apr 9, 2015 #3
    Yes my goal is to plot a potential vs. displacement curve. Is this even possible considering the number of functions in the first equation above?
    Here's some of my code: For simplicity, I assume period T(E) = 1

    %Define Functions

    %For a,b = 1

    %U = ax^2 + bx^3

    %Assuming T(E) = 1

    U = @(x) x^2+x^3

    fun = @(E,U) T/(sqrt(U-E))^-1

    int = @(E) integral(fun,0,U)

  5. Apr 9, 2015 #4


    User Avatar
    Gold Member

    I think something like this is a step closer, but I'm not convinced it's correct yet.

    Code (Text):

    U = @(x) a*x.^2 + b*x.^3;
    A = (2*pi*sqrt(2))^-1;
    I = [];

    for n =0:0.1:10
      fun = @(E) T./(sqrt(U(n)-E));
      I = [I, A*integral(fun,0,U(n))];

    plot(0:0.1:10, I)
    The difference here is that I loop through possible values for U and evaluate the integral once for each one.
  6. Apr 9, 2015 #5
    Thank you I'll try it. I'm wondering if it would be more beneficial/ shed more light if I plotted x as a function of time. In that case I solved lagrange
    L = K - U = .5mv^2 - x^2 - x^3 since U = x^2 + x^3
    the associated equation of motion is then ma = -3x^2 - 2x
    I would then need to solve this differential equation to plot x(t). However I put this equation into wolfram alpha, and got a massive expression out and it makes me thing that it might not be possible to model this response.

    Is this the case?
  7. Apr 10, 2015 #6


    User Avatar
    Gold Member

    You can solve the second order differential equation numerically in MATLAB. Since the ODE functions only work with first order equations, you'll have to rewrite this second order equation

    m \frac{d^2 x}{dt^2} = -3x^2-2x

    As an equivalent system of first order equations,

    m \dot{y_2} = -3y_1^2-2y_1
    y_2 = \dot{x}
    y_1 = x

    You'll also need to define the initial conditions of the problem.

    Take a look at the examples here:

    Save the following function file, which defines the equations above (assume m=1):
    Code (Text):

    function dy = anharmonic(t,y)
      dy = zeros(2,1);
      dy(1) = y(2);
      dy(2) = -3*y(1).^2 - 2*y(1);
    Then set the time span and initial conditions before calling ode45:
    Code (Text):

    tspan = [0 15];
    y0 = [0; 0.25];
    [t,y] = ode45(@anharmonic, tspan, y0);

    You might want to check my work, because the plot looks more harmonic than anharmonic to me.

    Attached Files:

    Last edited: Apr 10, 2015
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook