MATLAB Examples

# Matlab example with demo 'ab' in 07P mode

This demonstration uses the 'ab' example of AUTO to demonstrate the use of the Dynamical Systems Toolbox in the 07P (AUTO) mode. Input and output files are similar to the ones used in AUTO.

## Find the equilibrium of the system

The continuation run has to start from an equilibrium state, or a stable limit cycle condition. Most applications will start from an equilibrium condition. This equilibrium is also known as a quasi-steady state, or a trim condition. There are two ways to obtain such an equilibrium:

1. Run the simulation for a sufficient period of time, until the derivatives (of the states of interest) are equal to zero. This is the easiest option.
2. Use a trim routine to obtain a trim point.

Choose/guess initial conditions

```PAR0=[0,14,2;0,14,2]'; ts=[0,10]'; U0=[0,0]'; ```

Run simulation. The equations need to be defined in a separate function file abode.m to that of the continuation function file, so that we can use the ode45 solver.

```options=simset('MaxStep',0.1); [T,Y]=ode45(@(t,U) abode(t,U,ts,PAR0),ts,U0,options); ```

Plot the results out. We can see that the state derivatives are zero. This is a trivial example seeing that we could have checked this by hand, but the method is applicable for more complex models.

```h=plot(T,Y); set(h,'LineWidth',2); title('Time history for Demo: ab'); xlabel('Time (s)'); ylabel('U(1),U(2)'); snapnow; close(gcf); ```

## Adjust equations file to contain correct initial conditions

We can now use the final states of the simulation as starting conditions for the continuation runs. If the derivatives for these states are not zero, the continuation runs will fail. Update the stpnt.m function to reflect the initial conditions.

Create an auto object

```a=auto; ```

Define input and output files

```a.s.RunMode='07P'; a.s.FuncFileName='abmatfunc'; a.s.OutFileName='ab'; ```

## Compute stationary solutions

Generate an AUTO restart file by running AUTO.

Make sure no restart files are present

```delete s.* ```

Copy constants file

```copyfile('c.ab.1','c.ab'); ```

Run continuation

```runauto(a); ```
``` ---------------------- AUTO07 FOR MATLAB --------------------- USER NAME : ECOETZEE DATE : 26/10/2010 10:02:11 MODEL NAME : ab CONSTANTS FILE : c.ab < BR PT TY LAB PAR(01) L2-NORM U(01) U(02) 1 1 EP 1 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 1 33 LP 2 1.05739E-01 1.48439E+00 3.11023E-01 1.45144E+00 1 70 LP 3 8.89318E-02 3.28824E+00 6.88982E-01 3.21525E+00 1 90 HB 4 1.30900E-01 4.27187E+00 8.95080E-01 4.17704E+00 1 92 EP 5 1.51242E-01 4.36975E+00 9.15589E-01 4.27275E+00 Total Time 0.938E-01 > ```

Save data to output files - b.ab, d.ab, s.ab

```sv('matab1'); ```
```Saving fort.7 as b.matab1... done Saving fort.8 as s.matab1... done Saving fort.9 as d.matab1... done ```

## Plot stationary solutions

We have seen from the output to the MATLAB command window that there are two limit points and one Hopf bifurcation. We have written a small routine to extract the data from the fort.7 output file. You can try to write your own routine, or you can use the simple routine ploteq. Plot the results out to see what the bifurcation diagram looks like. We plot the continuation paramter PAR(2) against the L2-NORM parameter. Note the plotting conventions that are used. Solid lines are used for stable equilibrium states, and broken lines for unstable equilibrium states. The limit points (LP) indicate a qualitative change in the stability, while the Hopf bifurcation (HB) indicates a transition from a steady state to a limit cycle.

```ploteq('matab1','PAR(2)','L2'); snapnow; ```

## Compute periodic solutions

We want to do the continuation from the Hopf Bifurcation to determine the amplitude and frequency of the limit cycle as we change the parameter.

Copy the constants file

```copyfile('c.ab.2','c.ab'); ```

Run continuation

```runauto(a); ```
``` ---------------------- AUTO07 FOR MATLAB --------------------- USER NAME : ECOETZEE DATE : 26/10/2010 10:02:12 MODEL NAME : ab CONSTANTS FILE : c.ab RESTART FILE : s.ab < BR PT TY LAB PAR(01) L2-NORM MAX U(01) MAX U(02) PERIOD 4 30 6 1.21347E-01 4.08243E+00 9.83711E-01 6.30390E+00 2.29801E+00 4 60 7 1.18591E-01 3.76857E+00 9.97665E-01 8.25264E+00 3.66907E+00 4 90 8 1.14794E-01 3.11323E+00 9.99585E-01 9.97959E+00 6.25524E+00 4 120 9 1.06928E-01 2.51769E+00 9.99368E-01 9.62065E+00 9.53437E+00 4 150 10 1.05558E-01 2.14343E+00 9.99115E-01 9.33106E+00 1.60574E+01 4 180 11 1.05507E-01 1.91563E+00 9.99107E-01 9.29963E+00 2.89632E+01 4 210 12 1.05507E-01 1.81363E+00 9.99096E-01 9.29899E+00 4.39597E+01 4 240 13 1.05507E-01 1.76140E+00 9.99092E-01 9.29758E+00 5.89585E+01 4 270 14 1.05507E-01 1.72960E+00 9.99094E-01 9.29921E+00 7.39580E+01 4 300 EP 15 1.05507E-01 1.70820E+00 9.99093E-01 9.29914E+00 8.89577E+01 Total Time 0.376E+02 > ```

Save the data into different files

```sv('matab2'); ```
```Saving fort.7 as b.matab2... done Saving fort.8 as s.matab2... done Saving fort.9 as d.matab2... done ```

## Plot periodic solutions

Add the L2-NORM of the Hopf-bifurcation. This shows up as magenta in the plot below.

```ploteq('matab2','PAR(2)','L2'); snapnow; ```

## Two-parameter continuation, follow locus of lower limit point (Forward)

Trace out the locus of the lower limit point while stepping in a forward direction. We vary the first and the thrid parameters for the continuation, so restart from label 2 in this first continuation, hence we need to copy the restart file from the first run to the appropriate name.

```copyfile('s.matab1','s.ab'); ```

Copy the constants file

```copyfile('c.ab.3','c.ab'); ```

Run continuation

```runauto(a); ```
``` ---------------------- AUTO07 FOR MATLAB --------------------- USER NAME : ECOETZEE DATE : 26/10/2010 10:02:51 MODEL NAME : ab CONSTANTS FILE : c.ab RESTART FILE : s.ab < BR PT TY LAB PAR(01) L2-NORM U(01) U(02) PAR(03) 2 27 LP 6 1.35335E-01 2.06249E+00 5.00226E-01 2.00091E+00 2.50000E+00 2 100 EP 7 2.07738E-07 1.82583E+01 9.45157E-01 1.82339E+01 -2.74306E-01 Total Time 0.781E+00 > ```

Save the data

```sv('matab3'); ```
```Saving fort.7 as b.matab3... done Saving fort.8 as s.matab3... done Saving fort.9 as d.matab3... done ```

## Two-parameter continuation, follow locus of lower limit point (Backward)

Trace out the locus of the lower limit point while stepping in a backward direction. We vary the first and the thrid parameters for the continuation, so restart from label 2 in this first continuation, hence we need to copy the restart file from the first run to the appropriate name.

```copyfile('s.matab1','s.ab'); ```

Copy the constants file

```copyfile('c.ab.4','c.ab'); ```

Run continuation

```runauto(a); ```
``` ---------------------- AUTO07 FOR MATLAB --------------------- USER NAME : ECOETZEE DATE : 26/10/2010 10:02:52 MODEL NAME : ab CONSTANTS FILE : c.ab RESTART FILE : s.ab < BR PT TY LAB PAR(01) L2-NORM U(01) U(02) PAR(03) 2 35 EP 6 -1.31940E-03 9.96433E-01 -3.58652E-03 9.96426E-01 -1.05039E+00 Total Time 0.281E+00 > ```

Save the data

```sv('matab4'); ```
```Saving fort.7 as b.matab4... done Saving fort.8 as s.matab4... done Saving fort.9 as d.matab4... done ```

## Two-parameter continuation, follow locus of Hopf point (Backward)

Trace out the locus of the Hopf Bifurcation while stepping in a backward direction. We vary the first and the thrid parameters for the continuation, so restart from label 4 in this first continuation, hence we need to copy the restart file from the first run to the appropriate name.

```copyfile('s.matab1','s.ab'); % Copy the constants file copyfile('c.ab.5','c.ab'); ```

Run continuation

```runauto(a); ```
``` ---------------------- AUTO07 FOR MATLAB --------------------- USER NAME : ECOETZEE DATE : 26/10/2010 10:02:53 MODEL NAME : ab CONSTANTS FILE : c.ab RESTART FILE : s.ab < BR PT TY LAB PAR(01) L2-NORM U(01) U(02) PAR(03) 4 100 EP 6 8.81006E-05 1.17440E+01 9.14609E-01 1.17083E+01 9.36293E-02 Total Time 0.108E+01 > ```

Save the data

```sv('matab5'); ```
```Saving fort.7 as b.matab5... done Saving fort.8 as s.matab5... done Saving fort.9 as d.matab5... done ```

## Plot the locii of the limit points and Hopf bifurcation

Add the locii of the limit points and Hopf bifurcation

```ploteq('matab3','PAR(2)','L2'); ploteq('matab4','PAR(2)','L2'); ploteq('matab5','PAR(2)','L2'); snapnow; close(gcf); ```

## Plot the limit cycle response

We have written a small routine to extract the data from the fort.8 output file for the second run. Plot the response of the limit cycle behaviour at points 6, 8, 10. AUTO normalises the response so that the period is equal to 1 second. The plot on the left shows the response of state U(2), while the plot on the right shows a phase plot.

```plotlc('matab2'); snapnow; close(gcf); ```