This is a really complex example to start with if you're just getting into deriving equations of motion; I would think any course would have you do lots of examples with much simpler systems and build towards working on something like this.
That said, when I do work like this in Matlab I try to start with positional equations. You've got drawings there that show things like $c$, $h$, $\phi$, etc. Which of those are time-varying, and which aren't?
Perturb the system and make your positional equations from there, like the second drawing. Look at examples for compound pendulums and see how they've setup their positional equations. Remember that ultimately you're trying to go for kinetic and potential energies, so you'll need to get positions to all the centers of mass.
When you create the equations in Matlab, use the symbolic toolbox. You can define constants with the sym
command, like
sym b;
and time-varying parameters by putting (t)
after the variable, like
sym s(t);
then, when you write your equations, you can just use the variables b
and s
and Matlab understands that s
varies as a function of time and b
does not. This is important because, once you get your position equations completed, you can take derivatives with the diff
command. You can just do diff(myEqn)
and Matlab understands to take the first derivative with respect to t
, or you can be explicit about it like diff(myEqn, t)
. There's are a few other specifics you can read about in the documentation, but those are the basics for diff
.
The Euler-Lagrange equation has the "q" terms, which is just the "variable in question" at the time, so you repeat the process for each variable and wind up with your system of equations. I think you would use $s$ as a q once and then use $\phi$ as a q once, which should end you up with some set of equations for $\ddot{s}$ and $\ddot{\phi}$, which you'd use for your equations of motion. If you're going to put an actuator on $b$ or something to change the bob weight then you'd want to use that, too.
Ultimately you need to get the Lagrangian to start with, like
$$
L = K - U \\
$$
and you can do that with $U = mgh$ and $K = (1/2)mv^2$. You'll get your velocity $v$ term by taking the time derivative of your position equation, with diff
.
Then, when you get $L$, you can find the Euler-Lagrange equation, like
$$
\frac{\partial L}{\partial q} L(t) - \frac{d}{dt}\frac{\partial L}{\partial \dot{q}} L(t) = 0 \\
$$
So, once you define what $L$ is, you need to again go through each time for a different "q" variable ($s$, $\phi$, etc.) and get your Euler-Lagrange equation. Again, you can kind of leverage Matlab's diff
function for this but it becomes more difficult because Matlab does a pretty poor job at handling partial derivatives of derivatives.
What I do is to use the subs
command to pull my derivatives out and replace them with something easier to reference, like
syms sDot;
L2 = subs(L, diff(s,t), sDot);
That would give you L2
as a function of sDot
instead of diff(s, t)
. I do this because again the partial derivative is kind of a pain to do. What you're going for is taking the derivative of L
with respect to $\dot{s}$, but what that looks like in Matlab would be:
dLdqDot = diff(L, diff(s,t));
and it's been long enough that I can't remember exactly what the issue was, but I remember Matlab not being able to intelligently handle that. By swapping out your diff(s, t)
with a singular variable like sDot
you can instead do something like:
dLdqDot = diff(L, sDot);
and get better results. Once you're done with that, you'll need to subs
back in your version as a function of time so you can take the time derivative of that, the $d/dt$ of the $d/dt(\partial L / \partial \dot{q})$ term.
dLdqDot_t = subs(dLdqDot, sDot, diff(s, t));
Then you get your right hand side:
eulerLagrange_RHS = diff(dLdqDot_t, t);
and then finally one more round of subs
to get your variables back out to be sortable again:
syms sDdot;
el_RHS = subs(eulerLagrange_RHS, diff(s, t, t), sDdot);
el_RHS = subs(el_RHS, diff(s, t), sDot);
and it's important to note that you need to do the second derivative FIRST or you'll wind up with something like diff(sDot, t)
instead of diff(s, t, t)
or sDdot
. After you do the second derivative you can do the first derivative.
Then, once you're done with the (d/dt)(dL/dqDot)
then you can work on doing the same for the dL/dq
as the left-hand side, and then you can setup your equation:
equationOfMotion(1) = el_LHS - el_RHS == 0;
and then you can solve(equationOfMotion(1), sDdot)
and get your equation of motion for $\ddot{s}$. Then you do it all over again for $\phi$. And yes, it's super, super tedious, but again you should have had at least one full semester and probably two full semsters of doing homework basically like this at least once a week, classes like statics, dynamics, system dynamics, etc.
This is all more about how to wrangle Matlab than the specifics about your problem, but do your work in a script in Matlab and you can make all the pretty plots you want along the way to check your work.
I would strongly, strongly recommend you parameterize your variables at the end and use a variety of plot
commands to visualize what you're working with, to recreate the scene through various values of $\phi$, etc., to make sure you haven't gone wrong somewhere.
If you try working your way through this and have more issues, please come back and show your work with the problem as far as you've gotten and a description of why you think it's wrong or what has you stuck.
Also a disclaimer - I don't actually have access to a copy of Matlab any more or any of my prior work (changed jobs wooo) so this is all from memory; you'll probably find some mistakes in the syntax or something but the steps should be generally accurate.