Physics problem help!

We were given a problem in my physics class that we are supposed to code, despite not being given much on how to use matlab. I've been at this for hours and I can't seem to figure it out and any type of help will be appreciated.

Problem 3: Oscillation of a Sodium atom in a Sodium-Chloride crystal

Problem: The Sodium-Chloride crystal is an example of a ionic crystal. It consists of alternating positively and negatively charged ions, the Sodium and Chloride atoms respectively. A Sodium atom moves subject to forces from the surrounding ions - and an analysis of the system shows that for small perturbations from the equilibrium configuration, a Sodium ion will experience an acceleration of a = −kx , where k = 0.17×1027s^−2 , where x is the position of the ion measured in meters. Assume that the ion starts at position x0 = 0 time t = 0 with a velocity v0 = 364m/s. (a) Find the time development of x(t). (b) Find the maximum value of x. (c) Find the oscillation time, that is the time it takes until the ion is in the same state as it started in.

[James Tursa]

Supposed to code it how? Using your own written ordinary differential equation code, such as Euler method? Using a built-in integrator such as ode45? Symbolically in MATLAB? Solving by hand?

E.g., see this page for the close form solution:

https://en.wikipedia.org/wiki/Harmonic_oscillator

[Joe Worster]

Yes, such as the Euler method I believe

Answer 1

[James Tursa]

So, I will give you a start for an Euler scheme. You have two states to keep track of, the position and the velocity. Define two variables for this:

x is the position
v is the velocity

You are given two initial conditions for these states:

x0 = 0 m
v0 = 364 m/s

And a force constant

k = 0.17×10^27 s^−2   <-- I assume this is 10^27 ???

For an arbitrary variable Y, the Euler method of propagating is

Y(t+dt) = Y(t) + (dY/dt) * dt;

So putting the pieces together here is an outline of your code:

% Initial conditions
x0 = 0; % m
v0 = 364; % m/s
k = 0.17e27; % s^−2
n = 1000; % Number of steps to integrate (you will have to experiment with this value)
dt = ________; % You will have to experiment with this value also
% Allocate states
x = zeros(1,n);
v = zeros(1,n);
t = zeros(1,n);
% Initialize first state
x(1) = x0;
v(1) = v0;
t(1) = 0;
% Loop for the integration
for i=1:n-1
    x(i+1) = x(i) + ___________;
    v(i+1) = v(i) + ___________;
    t(i+1) = t(i) + ___________;
end

You need to figure out what goes in the blanks above. The x equation will be based on the derivative of x, which is velocity. The v equation will be based on the derivative of v, which is acceleration. The t equation is just to increment t by dt.

[Joe Worster]

Thanks for the help.

For my "n" value I just did n = ceil(time/dt);

For "dt" I ended up with dt = 1e-37 because I had to make it as small as possible due to my "time" variable being extremely small because we are doing the oscillation time of an atom. My time variable is time = 0.00000000000000000000000000000000000025; (professor told me to experiment with numbers and this is what finally allowed something to show up on my graph along with my dt variable).

In my for loop I ended up with this:

v(i+1) = v(i) + a(i)*dt;
x(i+1) = x(i) + v(i+1)*dt;
t(i+1) = t(i) + dt;

On the graph, it is just showing a straight, increasing line which is where I am starting to become confused because I believed this was supposed to be oscillating. Any suggestions? (Included picture of the graph)

[James Tursa]

You are not integrating enough time. Increase your time.

Also, typically one uses the current state for the derivatives in Euler's method. So the x equation would be:

x(i+1) = x(i) + v(i)*dt;

What are you using for a(i)?

[Joe Worster]

Whenever I increase my time variable my graph goes completely blank, or am I just not following very well?

Here's my current script:

k = 0.17e27; 
x0 = 0.0;  %initial position in meters
v0 = 364; %m/s
time = 0.00000000000000000000000000000000000025;
dt = 1e-37;
n = ceil(time/dt);
x = zeros(1,n);
v = zeros(1,n);
t = zeros(1,n);
x(1) = x0;
v(1) = v0;
t(1)=0;
for i = 1:n-1 %formulas
    a = -k*x;
    v(i+1) = v(i) + a(i)*dt;
    x(i+1) = x(i) + v(i)*dt;
    t(i+1) = t(i) + dt;
end

      plot(t,x)
      xlabel('t[s]');
      ylabel('v[m/s]');

[Joe Worster]

 if true
k = 0.17e27; 
x0 = 0.0;  %initial position in meters
v0 = 364; %m/s
time = 0.00000000000000000000000000000000000025;
dt = 1e-37;
n = ceil(time/dt);
x = zeros(1,n);
v = zeros(1,n);
t = zeros(1,n);
x(1) = x0;
v(1) = v0;
t(1)=0;
for i = 1:n-1 %formulas
    a = -k*x;
    v(i+1) = v(i) + a(i)*dt;
    x(i+1) = x(i) + v(i)*dt;
    t(i+1) = t(i) + dt;
end

    plot(t,x)
    xlabel('t[s]');
    ylabel('v[m/s]'); 
 end

figured I would repost that mess I just posted, just ignore the "if true" and "end" statements at the top and bottom of it.

[James Tursa]

It helps here to know the solution in advance to pick appropriate values for dt and n without a lot of trial and error. If you look at the wiki link I posted above, the period of a harmonic oscillator is:

period = T = 2*pi/sqrt(k/m)

In your case you are using an equivalent mass value of 1, so we just get

period = 2*pi/sqrt(k) = 2*pi/sqrt(0.17e27) = 4.82e-13

So, pick dt to make sure you have enough samples per period to make a meaningful integration and plot. This is somewhat arbitrary, but let's just make it 1000. Then

dt = period / 1000;

And if you want to run three full periods, then

n = 3000;

Give these values a try. Note that the period of 4.82e-13 is much larger than your time duration of 2.5e-37. Essentially you are currently only plotting a very small piece of the full period, and that very small piece looks linear.

P.S. Your y label should be 'x[m]', not 'v[m/s]'.

[James Tursa]

Another thing I noticed was these lines:

    a = -k*x;
    v(i+1) = v(i) + a(i)*dt;

This works, but it looks strange since the a = -k*x part is multiplying k by the entire x vector at each iteration. You don't really need this. All you need is to multiply k by the current x position. E.g.,

    a = -k*x(i);
    v(i+1) = v(i) + a*dt;

[Jim Riggs]

Did you notice that time is set to 2.5e-37 and dt = 1e-37.

Thus, the number of calculations, n = ceil(time/dt) = 3

The calculation lop is for i=1:n-1, so he is only doing two calculations.

[James Tursa]

@Jim Riggs: Yes, I noticed that. I first advised him to increase his total integration time hoping he would look at that himself, and the latest is I have actually given him appropriate integration times and step sizes to use. So, hopefully he can get everything fixed up now.

[Joe Worster]

Thanks so much for the help!