using MATLAB ODE solver for crane with hanging load

Question

I'm trying to numerically evaluate a coupled ODE using MATLAB. The system I'm moddeling is a crane (mass $m_t$) with a hanging load (mass $m_p$, Inertia $I$ at distance $l$ under the crane). I have a function $u(t)$ that applies a force in the $x$ direction (horizontal) on the crane.

I've already determined (and double checked) the ODE's, which are:

$$(I+m_pl^2)\ddot{\theta}+m_pgl\sin(\theta)=-m_pl\ddot{x}\cos(\theta)$$ $$(m_t+m_p)\ddot{x}+b\dot{x}+m_pl\ddot{\theta}\cos(\theta)-m_pl\dot{\theta}^2\sin(\theta)=u$$

Trying to put this in MATLAB I do this:

b = 0.1;
l = 1;
m_p = 1;
m_t = 1;
g = 9.81;
I = 1;

function f = u_tFn(t)
    f = zeros(size(t));
    f(t>0 & t<1) = 1;
end

function dstatedt = ode(t, state)
    theta = state(1);
    dthetadt = state(2);
    d2thetadt2 = (-m_p*g*l*sin(theta)-m_p*d2xdt2*l*cos(theta))/(m_p*l^2+I);
    x = state(3);
    dxdt = state(4);
    d2xdt2 = (-b*dxdt-m_p*l*d2thetadt2*cos(theta)+m_p*l*dthetadt^2*sin(theta)+u_tFn(t))/(m_t+m_p);
    dstatedt = [dthetadt;d2thetadt2;dxdt;d2xdt2];
end

But that doesn't work... As you can see, the problem is that $\ddot{\theta}$ is dependent on $\ddot{x}$ and vice versa $\ddot{x}$ is dependent on $\ddot{\theta}$.

How should I approach this problem?

Answer 1

Ok, so J. M. is not a mathematician helped me find an answer.

As he writes, you can seperate the system into linear equations. $$\begin{pmatrix}I+m_pl^2&m_pl\cos(\theta)\\m_pl\cos(\theta)&m_t+m_p\end{pmatrix}\begin{pmatrix}\ddot{\theta}\\\ddot{x}\end{pmatrix}=\begin{pmatrix}-m_pgl\sin(\theta)\\u+m_pl\dot{\theta}^2\sin(\theta)-b\dot{x}\end{pmatrix}$$

And of course $\begin{pmatrix}a&b\\c&d\end{pmatrix}^{-1}=\frac{1}{ab-cd}\begin{pmatrix}d&-b\\-c&a\end{pmatrix}$

So $$\begin{pmatrix}\ddot{\theta}\\\ddot{x}\end{pmatrix}= \frac{1}{I(m_p+m_t)+l^2m_p(m_p\sin^2(\theta)+m_t)} \begin{pmatrix}m_p+m_t&-lm_p\cos(\theta)\\-lm_p\cos(\theta)&m_pl^2+I\end{pmatrix} \begin{pmatrix}-m_pgl\sin(\theta)\\u+m_pl\dot{\theta}^2\sin(\theta)-b\dot{x}\end{pmatrix}$$

Which removes the cyclic dependencies. So I changed my MATLAB code to:

function dstatedt = ode(t, state)
    u = u_tFn(t);
    theta = state(1);
    dthetadt = state(2);
    x = state(3);
    dxdt = state(4);
    determinant = 1/(I*(m_p+m_t)+l^2*m_p*(m_p*sin(theta)^2+m_t));
    d2thetadt2 = determinant * (...
        -(m_p + m_t)*m_p*g*l*sin(theta)...
        -l*m_p*cos(theta)*(u-b*dxdt+m_p*l*dthetadt^2*sin(theta)));
    d2xdt2 = determinant * (...
        l*m_p*cos(theta)*m_p*g*l*sin(theta)...
        +(m_p*l^2+I)*(u-b*dxdt+m_p*l*dthetadt^2*sin(theta)));
    dstatedt = [dthetadt;d2thetadt2;dxdt;d2xdt2];
end

But then decided I can have MATLAB solve this for me (gives more readable code imho)

function dstatedt = ode(t, state)
    u = u_tFn(t);
    theta = state(1);
    dthetadt = state(2);
    dxdt = state(4);
    d2dt2 = [I+m_p*l^2 m_p*l*cos(theta);...
             m_p*l*cos(theta) m_t+m_p]\...
        [-m_p*g*l*sin(theta);u-b*dxdt+m_p*l*dthetadt^2*sin(theta)];
    dstatedt = [dthetadt;d2dt2(1);dxdt;d2dt2(2)];
end

Still, I very much appreciate the hints on how to approach this problem. Thank you!