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

: ,

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

## 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 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 is fixed.

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

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 . 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); 