Get documents about "MatCont Worksheet Case studies on the numerical analysis of parameter
Document Sample
MatCont, Worksheet 5: Case studies on the numerical analysis
of parameter dependent ODEs
This worksheet poses two open-ended exercises and provides almost complete Matlab
scripts for solving them, the last part being left open for completion.
14 2.5
15
12
10 2
10
8
max(x1)
period
2
H
x
5 1.5
6
4
0 1
2
0
BP H 10 25
5 20
15 0.5
0 5 10 15 20 25 30 35 x 0 14 16 18 20 22 24 26
ρ 1 ρ ρ
Figure 10: Bifurcation diagram of the Lorenz equation (left), the branch of periodic
solutions (middle) and the period as a function of the parameter (right).
7. The Lorenz system
Investigate the Lorenz system
˙
x = σ(y − x),
˙
y = ̺x − y − xz,
˙
z = xy − βz,
for parameter values in the range ̺ ∈ [0, 30], σ = 10 and β = 8/3.
The bifurcation diagram is shown in the left-hand panel of Figure 10. The Lorenz
system has an equilibrium point at (x, y, z) = (0, 0, 0) for all values of ̺. A pitchfork
bifurcation occurs on this branch of equilibria at ̺ = 1. Switching to the intersecting
branch of equilibria we find a Hopf-bifurcation point near ̺ ≈ 25. The periodic solu-
tions emerging from this Hopf-point are shown in the middle panel of Figure 10. They
approach a homoclinic orbit from the equilibrium (x, y, z) = (0, 0, 0) to itself as the
parameter ̺ is decreased. The period of these periodic solutions tends to infinity as they
come closer and closer to this homoclinic orbit as can be seen in the right-hand panel of
Figure 10.
The Matlab script used for computing the data shown in Figure 10 is printed below.
It is complete except for the last continuation of the limit cycles that was performed
with increased accuracy.
a. Repeat the computations step by step using the script below as a template.
b. Complete the script with the missing continuation starting at the last limit cycle
computed. Use the function init LC LC to initialise this continuation and set
ntst = 50 and ncol = 4. Try also other values for ntst and ncol. What is the
lowest value of ̺ and the highest value of the period that you can reach?
29
c. The famous Lorenz attractor appears for values of ̺ to the right of the Hopf-
bifurcation point. Start a long-time simulation (that is, tspan=[0 100] or longer)
near the origin (x, y, z) = (0, 0, 0) for ̺ = 27 and visualise the steady-state.
global cds lds;
clf;
dim = 3; % dimension of ODE
pars = [10; 0; 8/3]; % = [si ro be]
ap = 2; % active parameter ro = #2
ode = @lor;
% equilibrium
if exist(’ep1.mat’, ’file’)~=2
opt=contset;
opt=contset(opt,’MaxNumPoints’,100);
opt=contset(opt,’Singularities’,1);
[x0 v0]=init EP EP(ode, [0;0;0], pars, ap);
[x v s h f]=cont(@equilibrium,x0,v0,opt);
save(’ep1’);
else
load(’ep1’);
end
cpl(x,v,s,[dim+1 1]);
sbp = s(2); % start data BP->EP
idx = sbp.index;
xbp = x(1:end-1,idx);
pbp = pars;
pbp(ap) = x(end,idx);
% branch of equilibria BP->EP
if exist(’ep2.mat’, ’file’)~=2
opt=contset(opt,’MaxNumPoints’,500);
opt=contset(opt,’MaxStepsize’,0.1);
opt=contset(opt,’InitStepsize’,0.05);
[x0 v0]=init BP EP(ode, xbp, pbp, sbp, 0.01);
[x v s h f]=cont(@equilibrium,x0,v0,opt);
save(’ep2’);
else
load(’ep2’);
end
cpl(x,v,s,[dim+1 1]);
sh = s(2); % start data H->LC
idx = sh.index;
xh = x(1:end-1,idx);
30
ph = pars;
ph(ap) = x(end,idx);
% branch of limit cycles H->LC
if exist(’lc1.mat’, ’file’)~=2
opt=contset(opt,’MaxNumPoints’,30);
opt=contset(opt,’MaxStepsize’,5);
opt=contset(opt,’InitStepsize’,1);
[x0 v0]=init H LC(ode, xh, ph, ap, 0.01, 10, 4);
[x v s h f]=cont(@limitcycle,x0,v0,opt);
save(’lc1’);
else
load(’lc1’);
end
[m points] = size(x);
xx = x(1:end-2,:);
xx = reshape(xx, [dim (m-2)/dim points]);
xa = max(squeeze(xx(1,:,:)), [], 1);
xi = min(squeeze(xx(1,:,:)), [], 1);
p1 = x(end,:);
cpl([p1; xa],v,s,[1 2]);
cpl([p1; xi],v,s,[1 2]);
8. An inverted pendulum (a reduced system)
Investigate the equations of an inverted pendulum
˙
x = y,
˙
y = z,
z = −αx + γy + βz + x3 ,
˙
for parameter values in the range α = sin ϕ, β = cos ϕ cos ψ and γ = cos ϕ sin ψ, ϕ ∈
[0, π/2] and ψ = 4.4.
The bifurcation diagram is shown in the left-hand panel of Figure 11. The inverted
pendulum has an equilibrium point at (x, y, z) = (0, 0, 0) for all values of ϕ. A Hopf-
bifurcation occurs on this branch of equilibria near ϕ = 0.3. A branch of so-called
symmetric limit cycles emerges from this point and a pitchfork bifurcation (branch point
of limit cycles) occurs near ϕ = 0.5. A branch of so-called non-symmetric limit cycles
emanates as shown in the right-hand panel of Figure 11, the bifurcation breaks the
symmetry.
The Matlab script used for computing the data shown in Figure 11 is printed below.
It is complete except for the last continuation of the limit cycles that emerge in a
pitchfork bifurcation.
31
1
PD 1
0.5 BPC 0.5
max(x1)
2
0
x
H
0
−0.5
PD
−0.5
BPC −1
0.2 1
0.4 0
0.6
−1
0.8
−1 x
0 0.2 0.4 0.6 0.8 1 1.2
φ
1
φ
Figure 11: Bifurcation diagram of the inverted pendulum (left) and the branches of peri-
odic solutions (right).
a. Repeat the computations step by step using the script below as a template.
b. Complete the script with the missing continuation of limit cycles emerging at the
branch-point of limit cycles (BPC). Use the function init BPC LC to initialise this
continuation. What bifurcations are detected along this branch?
c. Compute this branch also backwards.
global cds lds;
clf;
dim = 3; % dimension of ODE
pars = [0.1; 4.4]; % = [phi psi]
ap = 1; % active parameter phi = #1
ode = @ipen;
% equilibrium
if exist(’ep1.mat’, ’file’)~=2
opt=contset;
opt=contset(opt,’MaxNumPoints’,20);
opt=contset(opt,’Singularities’,1);
[x0 v0]=init EP EP(ode, [0;0;0], pars, ap);
[x v s h f]=cont(@equilibrium,x0,v0,opt);
save(’ep1’);
else
load(’ep1’);
end
cpl(x,v,s,[dim+1 1]);
% branch of limit cycles H->LC
if exist(’lc1.mat’, ’file’)~=2
idx = s(2).index;
32
x1 = x(1:end-1,idx);
pars(ap) = x(end,idx);
opt=contset(opt,’MaxNumPoints’,100);
opt=contset(opt,’MaxStepsize’,0.1);
opt=contset(opt,’InitStepsize’,0.05);
[x0 v0]=init H LC(ode, x1, pars, ap, 0.001, 20, 4);
[x v s h f]=cont(@limitcycle,x0,v0,opt);
save(’lc1’);
else
load(’lc1’);
end
[m points] = size(x);
xx = x(1:end-2,:);
xx = reshape(xx, [dim (m-2)/dim points]);
xa = max(squeeze(xx(1,:,:)), [], 1);
xi = min(squeeze(xx(1,:,:)), [], 1);
p1 = x(end,:);
cpl([p1; xa],v,s,[1 2]);
cpl([p1; xi],v,s,[1 2]);
33
Related docs
