Biocybernetics lab
Polly Vacher 184 on Mon 22-11-2021
We are using the non-linear pendiulm equations from lab 1 to investigate ways of looking at non linear systems.
The equation for a simple pendulum with damping is
\[ \tau=ml\ddot\theta+mg\sin\theta+B\dot\theta \]where $m$ is the mass of the bob, $\theta$ is the angle between the gravity vector $g$ and the pendulum and $B$ is a linear term to represent friction and damping, and $\tau$ is the torque at the pivot. Mostly we will set $\tau=0$ and $B=0$.
Do a direct solve of the non-linear undamped ($B=0$) pendulum equations and plot the direct form to show the phase plane. For this you will need to create a range of values for the angle variable theta.
You can use radians or degrees. You will need to cover about 3 revolutions so either
>> theta=-10:.1:10 % theta measured in radians
>> theta=-700:700 % theta measured in degrees
You need to calculate the non-linear pendulum formula
\[ \omega= \pm\sqrt{ \frac{g}{l}\left(\cos\theta-\cos\theta_0\right) +\omega^2_0 } \]Most of these values are constants so set up
>> g=9.8;
>> l=somevalue
>> theta0=somevalue;
>> omega0=somevalue;
You can then calculate the function above and plot
>> plot(theta,omega)
do not forget that a square root function has both a positive and a negative answer!
Now do a numerical solve of the non-linear undamped pendulum equations. These equations are
\begin{align} \dot\theta&=\omega\\ \dot\omega&=-\frac{g}{l}\sin\theta \end{align}And like Quest 2, we need to set up some constants, so
>> g=9.8;l=somevalue
We also need theta0
and omega0
but this time as the starting value for the numerical integration.
>> theta0=somevalue; omega0=somevalue;
Finally the function to integrate. This function simply restating the above equation.
>> pendulumfn=@(t,x) [x(2); (g/l)*sin(x(1))];
You can test it with some typical values, for example
>> pendulumfn(42,[0;3])
Finally integrate with the starting values for some period of time.
>> [t,y]=ode45(pendulumfn,[0 20],[theta0;omega0]);
This time you will plot the angle (y(1)) against the angular velocity (y(2))
>> plot(y(:,1),y(:,2))
Try this plot with other starting values for theta0
and omega0
.
Compare your results with Quest 2
Now do a numerical solve of the non-linear damped pendulum equations. The equations above are simply augmented to
\begin{align} \dot\theta&=\omega\\ \dot\omega&=-\frac{g}{l}\sin\theta-\beta \omega \end{align}The only thing that needs to change is to update the pendulum function
>> pendulumfn=@(t,x) [x(2); (g/l)*sin(x(1))-beta*x(2)];
Try positive and negative values of beta and different values for omega0
and plot as before.
We will choose a simple model that is difficult to fit with least squares, that is
\[ y=A\sin(b x) \]where A and b are the model parameters and y and x are the set of data. These parameters will be considered as variables in the input vector $u$ so that$u=\begin{bmatrix}A&b\end{bmatrix}$ will be the input to the function we wish to minimise $err=f(u)$ s
So generate some data, to keep things simple keep x between 0 and 10 ($0\le x\le10$). Make the data look line a sin wave but you can use any number of cycles, any amplitude and any number of points. It is helpful to keep the numbers in sequence but only because that means that you can plot them with lines rather than as a scatter plot.
>> x=[....] % values between 0 and 10
>> y=[...] % any set of real number (positive or negative) should look `sinwaveish'
>> plot(x,y)
Now create the function $f(\vec{u})$ to minimise. $x$ and $y$ are now considered constants, and the parameters as variables. We can do this by creating a Matlab file with the editor and beginning the file with the word function
>> edit sinfittingfn.m
This function needs three lines, the first identifies the function, the second to calculate the error (for each pair $x,y)$ and the third to return the sum of error squared
Test the file by running with random numbers
>> sinfittingfn([1.1 3],x,y)
This should give you a single number that is the sum of squares error over the data for a particular model, in this case the model $y=1.1sin(3x)$
Once you can run the file without errors fit the function from a starting point.
>> [Usoln,fval] = fminsearch(@(u) myfun(u,x,y),[3 10]) %starting solution from A=3 and b=10
Now test
>> tx=0:.05:10;
>> plot(x,y,tx,Usoln(1)*sin(Usoln(2)*tx))
fval
?)Last modified 22/11/2021