Ok I have a project to do and I have to solve it numerically using
some approximation method in a program such as matlab, Maple, etc... 
My biggest stumbling block is just entering the differential equation
correctly. My question- how do you enter this equation and plot it
numerically: (xdot=sqrt((k/m)*x0^2-x^2)) xdot is the derivative of x,
x0 is just the initial x value and k and m are the spring constant and
mass in the system.  This equation is for a simple horizontally
oscillating mass srping system, the equation above is one my professor
wrote to describe it, and it should have a graph that looks like
x0*cos(x) when it is plotted.  how do you plot this numerically?

---

Request for Question Clarification by elmarto-ga on 19 Mar 2003 09:53 PST

Hello Mikek142!
I don't know about physics, but I'm teaching a differential equations
course for Economics. I know how to plot differential equations in
Matlab. However, I've tried plotting the function you asked for (with
arbitrary values for x(0), k and m), and it doesn't resemble
x(0)*cos(x). It looks something like log(x). I think it maybe has to
do with k and m, as I have no idea what a spring constant is and how
it relates to mass. So, is it possible that the graph doesn't look
like cos(x)? Can you tell me specific values for k, m and x(0).

Best regards,
elmarto-ga

---

Clarification of Question by mikek142-ga on 19 Mar 2003 16:38 PST

hey elmarto- im not sure if it was solved right.  try graphing the
first equation below (that's where it comes from) and see if it is a
sinusoid. the other equation i wrote should be similar

thanks a lot racecar- that matlab code certainly gives the right
answer. and it just shows me how absolutely little i know about
matlab- i would have propably not gotten that!!  and to clarify- my
professor wanted to write the equation of motion in a form where there
would be no second order term to make it easier to solve for
numerically, so he used the law of conservation of energy instead of
newton's second law of motion which gives you  m*x(doubledot)+k*x=0,
he took the total energy in the system at any instant to be equal to
the initial potential energy from the spring's initial displacement. 
this is true because it is a conservative system with no friction.  so
(1/2)k*x0^2=(1/2)k*x^2+(1/2)m*xdot^2
where the first term on the right is the potential energy of the
spring at any moment, and the second term on the right is the kinetic
energy of the block. (x0 is initial displacement) these must equal the
initial potential energy from the spring. that second equation should
be valid- he only wrote it that way to make it easier to put into a
numerical solver (actually it better be valid or it would violate the
laws of physics). since you told me how to solve the second order
system out right, I don't need it, so THANKS!!

(if you have the time, can you tell me what that xdot=zeros command is
about, and also in the command x0 = [1 0]'; , what does the ' do at
the end? and in the plot command, whats the colon for?

the actual equation I am trying to solve is
 M*g*L*(cos(Q)-cos(a))=(I/2)*(dQ/dt)^2  
which uses the same principle as in the spring example above except it
is a compound pendulum. it should also be oscillating. it is saying
the total energy in the system at any time must be equal the initial
potential engery in the pendulum by releasing it at alpha degrees
above equillibrium at time zero. Q is the angle of the pendulum,
theta, a=alpha is the initial angle of the pendulum from the
equillibrium (down) position, and I is the moment of inertia of the
pendulum. if anyone wants to check this out and see that this works in
giving an oscillating motion, go for it. ill slowly try to figure it
out too)

---

Clarification of Question by mikek142-ga on 29 Mar 2003 18:25 PST

alright I could not get it to work- how can you plot
dQ/dt=sqrt(((2*m*g*L)/I)*(cos(Q)-cos(a))) as described above.  I keep
getting "square root of a negative number" errors in Maple and Matlab,
or some weird wavy line that goes up and right. for a larger initial
angle, a, the period is supposed to become larger.

The equation isn't right.  For a mass m and a spring with spring
constant k, a correct equation is:

xdot = - sqrt(k/m) * x

provided that x = 0 is the equilibrium position.  (xdot means dx/dt).
If the initial position is x0 and the initial speed is zero, the
solution is

x = x0 cos(t * sqrt(k/m)).

Oops, I mean xdot means d^2x/dt^2.

And oops, get rid of the square root:

d^2x/dt^2 = -kx/m

That's the last oops.

To solve numerically in MATLAB, first define the following function:

%%%%%
function xdot = mikek(t,x)
xdot = zeros(2,1);
xdot(1) = x(2);        %xdot(1) = dx/dt
xdot(2) = -k/m*x(1);   %xdot(2) = d^2x/dt^2
%%%%%

Put in numbers for k and m.

Define a time interval and initial conditions, and call the function with ode45:

%%%%%
tspan = [0 8*pi/sqrt(k/m)];   % This choice will give you four periods
x0 = [1 0]';
[t,x] = ode45('mikek',tspan,x0);
plot(t,x(:,1));
%%%%%

The equation in the original question needs one more set of
parentheses:

xdot = sqrt((k/m)*(x0^2 - x^2))

To get the correct solution from this equation you have to choose xdot
positive half the time, and negative the other half.  The same is true
of your pendulum equation (I think it's for a simple pendulum, not a
compound one): because the equation only contains (dQ/dt)^2, you have
to guess the sign of dQ/dt.

The command xdot=zeros(2,1) defines xdot as a 2x1 column vector of
zeros.  The zeros don't matter, it's just there to set the size of the
vector.  In matlab, ' means transpose.  x0 needs to be a column
vector, not a row vector, and x0=[1 0] would give you a row vector.
The Nx2 matrix x has two columns: the first is the values of x
(displacment), and the second is xdot.  x(:,1) gives you the first
column.  It's a shorter way of writing x(1:end,1)--just means take all
the values in the first column.

im posting this again so the new comment will be shown on the first
page-

alright I could not get it to work- how can you plot
dQ/dt=sqrt(((2*m*g*L)/I)*(cos(Q)-cos(a))) as described above in the
clarification section.  I keep getting "square root of a negative
number" errors in Maple and Matlab,or some weird wavy line that goes
up and right. for a larger initial angle, a, the period is supposed to
become larger.