Two second order differential equation

[people]

hi, everybody
i got for homework a pendulum on a cart.

i solved the system and got two equations.

(m+M)x'' + M L theta'' = F

x'' + L theta'' + g* theta = 0

F = 1000 N, m = 500 kg, M = 1250 kg, l = 10m

i know how to put them in state space and solve them with SS block in
simulink.but i don't know how to solve the system with ode23. there are
no initial conditions stated in a task, i guess there are zero.

i understand this sintax
function xdot = vdpol(t,x)
xdot = [x(1).*(1-x(2).^2)-x(2); x(1)]

t0 = 0; tf = 20;
x0 = [0 0.25]'; % Initial conditions
[t,x] = ode23('vdpol',t0,tf,x0);
plot(t,x)

but don't know how to use it on those eq.

 

[matematikawan]

Do you know how to solve higher order differential equation using Runge-Kutta method?
Change your system of equations to the following form first before you use ode23

\frac{d\vec{y}}{dt} = \vec{f}(t,\vec{y})

where
\vec{y} = ( x , \dot{x} , \theta , \dot{\theta})^t

 

[people]

thank you for your reply.
i don't know if that is the method i am using but i know how to transform
this second order to first order eq.
and i got.
\ddot{\theta} +\frac{F}{m*l}+ \frac{(m+M)*g}{l*m}\theta=0
\ddot{x}-\frac{F}{m}-\frac{g*M}{m}\theta=0

x_{1} = \dot{x} -> \dot{x}_{1}= \ddot{x}
x_{2} = x -> \dot{x}_{2} = \dot{x} = x_{1}
\theta_{1} = \dot{\theta} -> \dot{\theta}_{1}= \ddot{\theta}
x_{2} = \theta -> \dot{\theta}_{2} = \dot{\theta} = \theta_{1}

and now i transform in
\dot{\theta}_{1} = -\frac{F}{m*l}-\frac{(m+M)*g}{l*m}\theta_{2}
\dot{\theta}_{2} = \theta_{1}
\dot{x}_{1} = \frac{F}{m}+\frac{g*M}{m}\theta_{2}
\dot{x}_{2} = x_{1}

and now i know how to put that in matlab
function thetadot = thetam(t, theta)
thetadot = [-0.2 - theta(2).*3.4335; theta(1)] % here are numbers because matlab
% can't "see" the variables F,m,M...
t0 = 0;
tf = 10;
theta0 = [0 0]';
[t,theta] = ode23('thetam',t0,tf,theta0);
plot(t,theta)
plot(t, theta(84:1:166))

and i get two graphs. one of them is the on i need.
can you explain me what is the second one, please.

but i don't know how to write in MATLAB solution for x
becase it has a theta in it eq??

 

[matematikawan]

try checking the size for theta. I think it's size is n by 2. That's why you have two graphs I think. It is just like using the command plot(t,theta(:,1),t,theta(:,2)).

The equation that you solve with MATLAB actually can be solve analytically because it is an inhomogeneous linear DE with constant coefficients.

 

[people]

yes, i know it has two columns. what is in a first? also solution to the eq?
problem is i don't know how to solve it for x... can you say exactly how
because i have no clue

 

[matematikawan]

My guess is that your second graph is t against theta. You must have known better. You wrote the program.

Try using the substitution that I wrote in my first post.

y_1= x , y_2=\dot{x},y_3=\theta , y_4=\dot{\theta}

 

[people]

i am sorry, but i can't see how is your's supstitution different form
mine when i wrote what i did to that point

 

[matematikawan]

Not much different actually . But you use subscripts for both x and theta that confuse me. It looks like you want to solve the two equations separately. What I have in mind is solving the two equations simultaneously using the four variables y₁, y₂, y₃, and y₄.

Then the plots for your solution are plot(t,y(:,1)) for x and plot(t,y(:,3)) for theta.

 

[people]

i forgot to write the solution... if anyone was interested or have similar problem
y_{1} = \dot{x}
y_{2} = x
y_{3} = \dot{\theta}
y_{4} = \theta
\dot{y}_{1} = \ddot{x}
\dot{y}_{2} = \dot{x} = y_{1}
\dot{y}_{3} = \ddot{\theta}
\dot{y}_{4} = \dot{\theta} = y_{3}

\dot{y}_{1} = \frac{F}{m} + \frac{gM}{m}y_{4}
\dot{y}_{3} = -\frac{F}{ml} - \frac{(m+M)g}{ml}y_{4}
\dot{y}_{4} = y_{3}
\dot{y}_{2} = y_{1}

and when we put it in matlab
_________________________
function ydot = ydot(t, y)
ydot = [2 + y(4).*24.525; y(1); -0.2 - y(4).*3.4335; y(3)]
%this is ydot.m file
_________________________
M = 1250;
m=500;
l = 10;
g = 9.81;
F = 1000;
t0 = 0;
tf = 10;
poc0 = [0 0 0 0]';
[t,rj]= ode23('xm', t0, tf, poc0)
plot(t,rj)

problem is similar to http://www.myphysicslab.com/pendulum_cart.html#indirect", except force is constant in my problem.