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

[baouba]

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?

Thanks

 

[kreil]

There is some useful info here:
http://arxiv.org/pdf/math-ph/0409034.pdf

I take it your goal here is to plot a potential vs. displacement curve? Any other info you can provide?

 

[baouba]

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)

ezplot('(-x+(sqrt(2)*2*pi)^-1)*int')

 

[kreil]

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

a=1;
b=1;
T=1;
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))];
end

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.

 

[baouba]

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?

 

[kreil]

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:
http://www.mathworks.com/help/matlab/ref/ode45.html

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

function dy = anharmonic(t,y)
  dy = zeros(2,1);
  dy(1) = y(2);
  dy(2) = -3*y(1).^2 - 2*y(1);
end

Then set the time span and initial conditions before calling ode45:

tspan = [0 15];
y0 = [0; 0.25];
[t,y] = ode45(@anharmonic, tspan, y0);
plot(t,y(:,1),t,y(:,2))
legend('Position','Velocity')

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