Biocybernetics lab

Polly Vacher 184 on Mon 22-11-2021

1. Non-linear pendulum equations

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$.

Part A: Direct solve

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!

Part B: Numeric solution of an undamped pendulum

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

Part C: Numeric solution of a damped pendulum

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.

2. Non-linear model fitting

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))
  1. Try different starting points, do you get the same model.
  2. Try using more or fewer data points.
  3. Note the effects of different local minima (what error is given by fval ?)

Last modified 22/11/2021