MATLAB Examples

Stationary States (1D Problem)(Demo : pd1)

This demo uses Euler's method to locate a stationary solution of a nonlinear parabolic PDE, followed by continuation of this stationary state in a free problem parameter. The equation is

: $\frac{\partial{u}}{\partial{t}}=D\frac{\partial^2{u}}{\partial{x^2}}+p_1u(1-u)$,

  • on the space interval [0,L] ,
  • where L=PAR(11)=10 is fixed throughout,
  • as is the diffusion constant D=PAR(15)=0.1.
  • the boundary conditions are u(0) = u(L) = 0 for all time.

Euler time integration is only first order accurate, so that the time step must be sufficiently small to ensure correct results. Indeed, this option has been added only as a convenience, and should generally be used only to locate stationary states.

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 pd1
%
f=[];
o=[];
dfdu=[];
dfdp=[];

f(1)= par(1) * u(1) * ( 1. - u(1) );


Definition of initial conditions file

Display initial conditions file contents.

type(a{1}.s.StpntFileName);
function [par,u,o]= stpnt(t)
%
% starting point for demo int
par=zeros(36,1);
o=[];

% set the (constant) parameter
par(1) = 1.;

% set the actual width of the space interval [0,par(11)]
par(11) = 10.;

% set the initial data in the (scaled) interval [0,1]
u(1) = sin(pi*t);

% also set the space derivative of the initial data
% note the scaling by 1/par(11) !
u(2) = pi * cos(pi*t)/par(11);

% set the diffusion constant
par(15) = 0.1;


Set intial conditions

In this case we load data from the starting point file. The stpnt.m file is called repeatedly for this option of Ips

Note that in the subroutine stpnt.m the initial data must be scaled to the unit interval, and that the scaled derivative must also be provided; see the equations-file func.f.

Initial data are $u(x)=\sin(\pi x/L)$ at time zero.

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

Load and display constants

In the first run the continuation parameter is the independent time variable, namely PAR(14), while $p_1=1$ is fixed.

The constants DS, DSMIN, and DSMAX then control the step size in space-time, here consisting of PAR(14) and $u(x)$.

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

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

  autoconstants handle

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


Time integration towards stationary state

Find stationary state.

a{1}=runauto(a{1});
Warning: Truncating U0 to length defined by Ndim 
 
    --------------- DYNAMICAL SYSTEMS TOOLBOX ---------------------     
 
USER NAME      : ECOETZEE
DATE           : 26/10/2010 10:10:47
 
 
<
  BR    PT  TY  LAB     TIME      L2-NORM     MAX U(01)
   1     1  EP    1   0.00000E+00   7.07107E-01   1.00000E+00
   1    50        2   2.29809E+00   8.80106E-01   9.97553E-01
   1   100        3   4.79305E+00   9.32012E-01   9.99496E-01
   1   150        4   7.29294E+00   9.38955E-01   9.99948E-01
   1   200        5   9.79295E+00   9.39752E-01   1.00005E+00
   1   205  EP    6   1.00430E+01   9.39773E-01   1.00007E+00

 Total Time    0.197E+01
>

Continuation of stationary states.

In the second run the continuation parameter is $p_1$. Restart from autof8 object from previous run

a{2}=a{1};
a{2}.c=cpd12(a{2}.c);
a{2}=runauto(a{2});
 
    --------------- DYNAMICAL SYSTEMS TOOLBOX ---------------------     
 
USER NAME      : ECOETZEE
DATE           : 26/10/2010 10:10:50
 
 
<
  BR    PT  TY  LAB      PAR(01)      L2-NORM     MAX U(01)
   1    25        7   6.45158E+00   9.76756E-01   1.00001E+00
   1    50        8   1.89456E+01   9.86502E-01   1.00000E+00
   1    75        9   3.14437E+01   9.89539E-01   1.00002E+00
   1   100  EP   10   4.39427E+01   9.91158E-01   1.00001E+00

 Total Time    0.252E+01
>

Plot the solution

Create plaut object and plot solution.

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