A Cusp Catastrophe with Sample Trajectories in SAS

Below is computer code written in SAS that graphs a cusp catastrophe surface, together with sample trajectories that move toward equilibria on that surface. Heun's method (a second-order Runge-Kutta) is illustrated for producing phase space trajectories, which is useful for introducing students to fourth-order Runge-Kutta methods.

options nocenter;

PROC IML;
* PARAMETER MEANINGS....
* 'T' is the level of social and political tension in the society...
it can be conceptualized as the proportion of the population
that is unhappy with the current state.
* 'a' is the sensitivity of partisan fragmentation to social Tension
* 'b' measures the response of change in fragmentation to existent
group dissonance (i.e., discontents meeting and grouping
together to form democratic partisan fragments).
* 'q' the sensitivity of partisan fragmentation to current levels of
partisan institutionalization (high levels of inst. should
slow down increasing fragmentation).
* 'L' this is the Lower bound for partisan institutionalization
* 'U' this is the Upper bound for partisan institutionalization
* 'p' is the sensitivity of change in institutionalization to
the relative level of Tension in society.
* 'e' measures the growth or loss rate of institutionalization based
on current levels of social Tensions and institutionalization.;

a= 1;b=-1;q= 1; T=0.5;
e=.5;L=-.8;U=.5;p=2.2;
h=.1;* h is the step size for Heun's method;

* NEWTON'S METHOD FOR THE EQUILIBRIUM ROOTS;
count1=0;
start;
do x1=-1 to 1 by 0.005;
do x0=-1 to 1 by 0.2;
x2=x0;
do j=1 to 5 by 1;
count1=count1+1;
deriv=(3#b#(x2##2)) + (a#T);
xnext=x2 - (((b#(x2##3)) + (a#T#x2) - (q#x1))/deriv);
xsave=xnext;xlast=x2;x2=xnext;
end;
* THE FILTER FOR THE ROOTS;
test=((b#(xnext##3)) + (a#T#xnext) - (q#x1));
abstest=abs(test);if abstest > .01 then goto skip1;
* END OF FILTER;
zero=(xnext||x1);zeros=zeros//zero;
skip1:
end;
end;

frag=zeros(|,1|);
inst=zeros(|,2|);

* HEUN'S METHOD FOR SAMPLE TRAJECTORIES;
do x4=-.9 to .9 by 0.6;x1=x4;
do x5=-1.5 to 1.5 by 0.5;x2=x5;x1=x4;
do k=1 to 50;
phmod2= ((b#(x2##3)) + (a#T#x2) - (q#x1));
phmod1= (-e#(T-(p#x1))#(x1-L)#(U-x1));
m2=x2+(h#phmod2);
m1=x1+(h#phmod1);
phmod2b= ((b#(m2##3)) + (a#T#m2) - (q#m1));
phmod1b= (-e#(T-(p#m1))#(m1-L)#(U-m1));
mod2next=x2+((h/2)#(phmod2 + phmod2b));
mod1next=x1+((h/2)#(phmod1 + phmod1b));
x2=mod2next;x1=mod1next;
F=F//x2;I=I//x1;end;
end;end;


Finish;Run;

allinst=I//inst;allfrag=F//frag;
trajects=allinst||allfrag;
party={'Inst' 'Frag'};
Create traject From trajects (|Colname=party|);
Append From trajects;
Close traject;
Run;
Data traject;Set traject;
Label inst='Control parameter A';
Label frag='X';
sym=1;

Symbol1 color=black v='.' f=simplex;
Symbol2 color=black f=simplex v='*';
Symbol3 color=black f=simplex v='.';
Symbol4 color=black f=simplex v='.';
Symbol5 color=black f=simplex v='.';
*TITLE1 f=simplex h=1.6 color=black
'Figure 1: A Singularity Model with';
*TITLE3 f=simplex h=1.6 color=black
'Partisan Fragmentation and Institutionalization';
title;

Proc Gplot Data=traject;
axis1 color=black minor=none
value=(h=1.5 f=simplex c=black)
label=(h=2 f=simplex c=black);
axis2 color=black minor=none
value=(h=1.5 f=simplex c=black)
label=(h=2 f=simplex c=black);
plot frag*inst=Sym/ nolegend skipmiss
vaxis=axis1 haxis=axis2 vminor=0 hminor=0;