Code covered by the BSD License  

Highlights from
Dynamical Systems Toolbox

image thumbnail

Dynamical Systems Toolbox

by

 

Bifurcation analysis of dynamical systems. Integration of AUTO bifurcation software into MATLAB.

Bratu's Equation (Demo : exp)

Bratu's Equation (Demo : exp)

This demo illustrates the computation of a solution family to the boundary value problem

: $u_1' = u_2$,

: $u_2' = -p_1  e^{u_1}$,

with boundary conditions $u_1(0)=0$, $u_1(1)=0$ This equation is also considered by Doedel, Keller, and Kernevez (1991).

Contents

Initialise workspace

Clear workspace

clear all

% Create a continuation object.
a{1}=auto;

Definition of function file

Display function file contents.

type(a{1}.s.FuncFileName);
function [f,o,dfdu,dfdp]= func(par,u,ijac)
%
% equations file for demo exp
%
f=[];
o=[];
dfdu=[];
dfdp=[];

f(1)= u(2);
f(2)=-par(1) * exp(u(1));


Definition of boundary conditions file

Display boundary conditions file contents.

type(a{1}.s.BcndFileName);
function [fb,o,dbc]= bcnd(par,u0,u1,ijac)
%
% boundary conditions for demo exp
%
fb=[];
o=[];
dbc=[];

fb(1)=u0(1);
fb(2)=u1(1);


Set intial conditions

We can either load data from the starting point file, or we can define the initial conditions directly into variable.

[a{1}.s.Par0,a{1}.s.U0,a{1}.s.Out0]=stpnt;

Load and display constants

Load the constants file.

a{1}.c=cexp1(a{1}.c);

% Display the constants.
a{1}.c
ans = 

  autoconstants handle

  Properties:
     Ndim: 2
    Noutx: 0
      Ips: 4
      Irs: 0
      Ilp: 1
      Icp: 1
     Ntst: 10
     Ncol: 4
      Iad: 3
      Isp: 1
      Isw: 1
     Iplt: 3
      Nbc: 2
     Nint: 0
      Nmx: 100
      Rl0: 0
      Rl1: 4
       A0: 0
       A1: 50
      Npr: 100
     Mxbf: 10
      Iid: 2
     Itmx: 8
     Itnw: 5
     Nwtn: 3
      Jac: 0
     Epsl: 1.0000e-007
     Epsu: 1.0000e-007
     Epss: 1.0000e-004
       Ds: 0.0100
    Dsmin: 1.0000e-003
    Dsmax: 0.2000
     Iads: 1
      Thl: []
      Thu: []
      Uzr: [6x2 double]


Compute the solution family containing the fold

Run equilbrium solutions.

a{1}=runauto(a{1});
 
    --------------- DYNAMICAL SYSTEMS TOOLBOX ---------------------     
 
USER NAME      : ECOETZEE
DATE           : 26/10/2010 10:10:02
 
 
<
  BR    PT  TY  LAB      PAR(01) INTEGRAL U(01)     MAX U(01)     MAX U(02)
   1     1  EP    1   0.00000E+00   0.00000E+00   0.00000E+00   0.00000E+00
   1    10  UZ    2   5.00071E-01   4.39344E-02   6.60050E-02   2.61316E-01
   1    13  UZ    3   1.00000E+00   9.32569E-02   1.40484E-01   5.49353E-01
   1    16  UZ    4   1.50000E+00   1.49861E-01   2.26457E-01   8.73235E-01
   1    19  UZ    5   2.00000E+00   2.16936E-01   3.28922E-01   1.24822E+00
   1    22  UZ    6   2.50000E+00   3.00803E-01   4.57926E-01   1.70436E+00
   1    26  UZ    7   3.00001E+00   4.17958E-01   6.39959E-01   2.31962E+00
   1    33  LP    8   3.51383E+00   7.61901E-01   1.18679E+00   4.00000E+00
   1    43  UZ    9   3.00000E+00   1.23967E+00   1.97495E+00   6.10338E+00
   1    49  UZ   10   2.49999E+00   1.50867E+00   2.43170E+00   7.20931E+00
   1    55  UZ   11   2.00000E+00   1.77497E+00   2.89390E+00   8.26876E+00
   1    62  UZ   12   1.50000E+00   2.07208E+00   3.41999E+00   9.42385E+00
   1    70  UZ   13   1.00000E+00   2.44387E+00   4.09144E+00   1.08469E+01
   1    82  UZ   14   5.00000E-01   3.01064E+00   5.13556E+00   1.29998E+01
   1   100  EP   15   1.48905E-01   3.89420E+00   6.80010E+00   1.63687E+01

 Total Time    0.286E+01
>

Restart at a labeled solution, using increased accuracy

a{2}=a{1};
a{2}.c=cexp2(a{2}.c);
a{2}=runauto(a{2});
 
    --------------- DYNAMICAL SYSTEMS TOOLBOX ---------------------     
 
USER NAME      : ECOETZEE
DATE           : 26/10/2010 10:10:06
 
 
<
  BR    PT  TY  LAB      PAR(01) INTEGRAL U(01)     MAX U(01)     MAX U(02)
   1   100  EP   16   9.52916E-05   8.40578E+00   1.56186E+01   3.40101E+01

 Total Time    0.364E+01
>

Plot the solution

Create plaut object and plot solution.

p=plautobj;
set(p,'xLab','Par','yLab','L2norm');
ploteq(p,a);

Contact us