Generate system files for Homoclinic RG flows¶
This Jupyter Notebook generates the system files for the model considered in the HomoclinicRGflows demo.
Add MatCont path and load sym package if GNU Octave is used¶
matcontpath = '../';
addpath(matcontpath);
addpath([matcontpath, '/Utilities']);
if isOctave
pkg load symbolic % for GNU Octave
end
Set the system name¶
system_name = 'HomoclinicRGflows';
Create coordinates and parameter names as strings¶
coordsnames = {'g1', 'g2', 'g3', 'g4'};
parnames={'epsilon','M', 'N'};
Create symbols for coordinates and parameters¶
The array par
is the array of symbols in the same order as parnames.
Due to the following two lines we may, for example, use either epsilon
or
par(1)
. There should no changes be need of this code.
syms(parnames{:}); % create symbol for alpha and delta
par=cell2sym(parnames); % now alpha1 is par(1) etc
syms(coordsnames{:}); % create symbol for alpha and delta
coords=cell2sym(coordsnames); % create 1 x n vector for coordinates
Define beta functions¶
beta12 = @(M,N) (1/8).*N.^(-2).*(96.*g2.^2.*g4.*((-8)+N).*N+768.*g2.*g3.*g4.*N.^2+ ...
128.*g1.*g3.*g4.*N.^2.*(3+5.*M+N+N.^2)+16.*g2.*g3.^2.*N.^2.*(16+ ...
7.*M+N+N.^2)+32.*g1.^2.*g4.*N.*((-40)+(-8).*M+8.*N+7.*M.*N+8.* ...
N.^2)+64.*g1.*g4.^2.*N.^2.*(22+(-2).*M+M.*N+M.*N.^2)+16.*g1.* ...
g3.^2.*N.^2.*(32+2.*M+N+M.*N+N.^2+M.*N.^2)+64.*g1.*g2.*g4.*N.*(( ...
-32)+(-4).*M+16.*N+2.*M.*N+5.*N.^2+2.*M.*N.^2)+4.*g2.^2.*g3.*N.*(( ...
-256)+(-64).*M+72.*N+10.*M.*N+19.*N.^2+2.*M.*N.^2)+16.*g1.^2.*g3.* ...
N.*((-80)+(-24).*M+30.*N+4.*M.*N+7.*N.^2+4.*M.*N.^2)+16.*g1.*g2.* ...
g3.*N.*((-144)+(-40).*M+42.*N+17.*M.*N+21.*N.^2+5.*M.*N.^2)+ ...
g2.^3.*(896+128.*M+(-352).*N+(-48).*M.*N+(-12).*N.^2+(-12).*M.* ...
N.^2+7.*N.^3+2.*M.*N.^3)+4.*g1.^2.*g2.*(928+224.*M+(-392).*N+( ...
-100).*M.*N+(-2).*N.^2+11.*M.*N.^2+23.*N.^3+7.*M.*N.^3+4.*N.^4)+ ...
4.*g1.^3.*(352+96.*M+(-152).*N+(-44).*M.*N+10.*N.^3+3.*M.*N.^3+M.* ...
N.^4)+2.*g1.*g2.^2.*(1600+320.*M+(-656).*N+(-136).*M.*N+32.*N.^2+ ...
10.*M.*N.^2+50.*N.^3+10.*M.*N.^3+5.*N.^4+2.*M.*N.^4)).*pi.^(-2);
beta22 = @(M,N) (1/8).*N.^(-2).*(1536.*g1.*g3.*g4.*N.^2+128.*g2.*g3.*g4.*N.^2.*(9+ ...
5.*M+N+N.^2)+32.*g1.*g3.^2.*N.^2.*(16+7.*M+N+N.^2)+192.*g1.^2.* ...
g4.*N.*((-12)+(-2).*M+4.*N+N.^2)+64.*g1.*g2.*g4.*N.*((-80)+(-16).* ...
M+16.*N+5.*M.*N+7.*N.^2)+64.*g2.*g4.^2.*N.^2.*(22+(-2).*M+M.*N+M.* ...
N.^2)+16.*g2.*g3.^2.*N.^2.*(48+9.*M+2.*N+M.*N+2.*N.^2+M.*N.^2)+ ...
16.*g1.^2.*g3.*N.*((-128)+(-32).*M+32.*N+12.*M.*N+14.*N.^2+3.*M.* ...
N.^2)+16.*g1.*g2.*g3.*N.*((-272)+(-72).*M+82.*N+15.*M.*N+24.*N.^2+ ...
6.*M.*N.^2)+16.*g2.^2.*g4.*N.*((-176)+(-40).*M+58.*N+14.*M.*N+17.* ...
N.^2+11.*M.*N.^2)+4.*g2.^2.*g3.*N.*((-576)+(-160).*M+176.*N+54.* ...
M.*N+74.*N.^2+17.*M.*N.^2)+8.*g1.^3.*(288+64.*M+(-112).*N+(-24).* ...
M.*N+(-6).*N.^2+3.*M.*N.^2+5.*N.^3+2.*M.*N.^3+N.^4)+2.*g1.*g2.^2.* ...
(3968+1024.*M+(-1600).*N+(-416).*M.*N+(-56).*N.^2+(-22).*M.*N.^2+ ...
85.*N.^3+17.*M.*N.^3+11.*N.^4)+4.*g1.^2.*g2.*(1856+448.*M+(-736).* ...
N+(-176).*M.*N+(-22).*N.^2+(-5).*M.*N.^2+43.*N.^3+8.*M.*N.^3+3.* ...
N.^4+2.*M.*N.^4)+g2.^3.*(2816+768.*M+(-1152).*N+(-320).*M.*N+12.* ...
N.^2+(-12).*M.*N.^2+69.*N.^3+30.*M.*N.^3+7.*N.^4+5.*M.*N.^4)).* ...
pi.^(-2);
beta32 = @(M,N) (1/8).*N.^(-3).*(384.*g1.*g2.*g4.*N.*(4+N.^2)+96.*g2.^2.*g4.*N.*( ...
8+N.^2)+128.*g1.*g3.*g4.*N.^2.*((-10)+(-2).*M+5.*N+M.*N+5.*N.^2)+ ...
32.*g3.^2.*g4.*N.^3.*(18+14.*M+7.*N+7.*N.^2)+64.*g3.*g4.^2.*N.^3.* ...
(22+(-2).*M+M.*N+M.*N.^2)+96.*g1.^2.*g4.*N.*(8+2.*N.^2+M.*N.^2)+ ...
24.*g3.^3.*N.^3.*(16+2.*M+2.*N+M.*N+2.*N.^2+M.*N.^2)+64.*g2.*g3.* ...
g4.*N.^2.*((-20)+(-4).*M+10.*N+2.*M.*N+5.*N.^2+2.*M.*N.^2)+16.* ...
g1.*g3.^2.*N.^2.*((-52)+(-14).*M+26.*N+7.*M.*N+14.*N.^2+7.*M.* ...
N.^2)+8.*g2.*g3.^2.*N.^2.*((-104)+(-28).*M+52.*N+14.*M.*N+38.* ...
N.^2+7.*M.*N.^2)+16.*g1.*g2.*g3.*N.*(208+56.*M+(-48).*N+(-12).*M.* ...
N+12.*N.^2+3.*M.*N.^2+8.*N.^3+M.*N.^3+2.*N.^4)+8.*g1.^2.*g3.*N.*( ...
208+56.*M+(-48).*N+(-12).*M.*N+12.*N.^2+6.*M.*N.^2+4.*N.^3+3.*M.* ...
N.^3+M.*N.^4)+8.*g1.^3.*((-96)+(-16).*M+24.*N+4.*M.*N+(-20).*N.^2+ ...
(-10).*M.*N.^2+6.*N.^3+2.*M.*N.^3+N.^4+M.*N.^4)+4.*g2.^2.*g3.*N.*( ...
416+112.*M+(-96).*N+(-24).*M.*N+18.*N.^2+6.*M.*N.^2+12.*N.^3+4.* ...
M.*N.^3+2.*N.^4+M.*N.^4)+g2.^3.*((-768)+(-128).*M+192.*N+32.*M.*N+ ...
(-96).*N.^2+32.*N.^3+4.*M.*N.^3+7.*N.^4+2.*M.*N.^4)+2.*g1.*g2.^2.* ...
((-1152)+(-192).*M+288.*N+48.*M.*N+(-176).*N.^2+(-40).*M.*N.^2+ ...
70.*N.^3+10.*M.*N.^3+21.*N.^4+2.*M.*N.^4)+4.*g1.^2.*g2.*((-576)+( ...
-96).*M+144.*N+24.*M.*N+(-104).*N.^2+(-40).*M.*N.^2+34.*N.^3+13.* ...
M.*N.^3+11.*N.^4+3.*M.*N.^4)).*pi.^(-2);
beta42 = @(M,N) (1/8).*N.^(-3).*(224.*g3.*g4.^2.*N.^3.*(2.*M+N+N.^2)+24.*g3.^3.* ...
N.^3.*(4+2.*M+N+N.^2)+224.*g1.*g4.^2.*N.^2.*((-4)+(-2).*M+2.*N+M.* ...
N+2.*N.^2)+16.*g1.*g3.^2.*N.^2.*((-16)+(-2).*M+8.*N+M.*N+7.*N.^2)+ ...
96.*g4.^3.*N.^3.*(8+(-2).*M+M.*N+M.*N.^2)+32.*g3.^2.*g4.*N.^3.*( ...
22+N+M.*N+N.^2+M.*N.^2)+224.*g2.*g4.^2.*N.^2.*((-4)+(-2).*M+2.*N+ ...
M.*N+N.^2+M.*N.^2)+64.*g2.*g3.*g4.*N.^2.*((-14)+(-4).*M+7.*N+2.* ...
M.*N+6.*N.^2+M.*N.^2)+64.*g1.*g3.*g4.*N.^2.*((-14)+(-4).*M+7.*N+ ...
2.*M.*N+2.*N.^2+2.*M.*N.^2)+8.*g2.*g3.^2.*N.^2.*((-32)+(-4).*M+ ...
16.*N+2.*M.*N+9.*N.^2+2.*M.*N.^2)+16.*g1.*g2.*g3.*N.*(80+16.*M+( ...
-12).*N+7.*N.^2+2.*M.*N.^2+N.^3)+8.*g1.^2.*g3.*N.*(80+16.*M+(-12) ...
.*N+4.*N.^2+2.*M.*N.^2+3.*N.^3+N.^4)+4.*g2.^2.*g3.*N.*(160+32.*M+( ...
-24).*N+17.*N.^2+7.*M.*N.^2+4.*N.^3+N.^4)+32.*g1.*g2.*g4.*N.*(96+ ...
36.*M+(-24).*N+(-12).*M.*N+2.*N.^2+(-2).*M.*N.^2+4.*N.^3+M.*N.^3+ ...
N.^4)+4.*g1.^3.*((-160)+(-48).*M+40.*N+12.*M.*N+4.*N.^2+2.*M.* ...
N.^2+4.*N.^3+M.*N.^3+2.*N.^4)+2.*g1.*g2.^2.*((-960)+(-288).*M+ ...
240.*N+72.*M.*N+(-88).*N.^2+(-20).*M.*N.^2+39.*N.^3+7.*M.*N.^3+ ...
13.*N.^4)+16.*g1.^2.*g4.*N.*(96+36.*M+(-24).*N+(-12).*M.*N+(-4).* ...
N.^2+(-2).*M.*N.^2+2.*N.^3+3.*M.*N.^3+M.*N.^4)+8.*g1.^2.*g2.*(( ...
-240)+(-72).*M+60.*N+18.*M.*N+(-8).*N.^2+(-1).*M.*N.^2+6.*N.^3+2.* ...
M.*N.^3+N.^4+M.*N.^4)+8.*g2.^2.*g4.*N.*(192+72.*M+(-48).*N+(-24).* ...
M.*N+10.*N.^2+2.*M.*N.^2+6.*N.^3+4.*M.*N.^3+N.^4+M.*N.^4)+g2.^3.*( ...
(-640)+(-192).*M+160.*N+48.*M.*N+(-96).*N.^2+(-24).*M.*N.^2+36.* ...
N.^3+12.*M.*N.^3+10.*N.^4+5.*M.*N.^4)).*pi.^(-2);
Define the system¶
dg1_dt = -epsilon*g1 + beta12(M,N);
dg2_dt = -epsilon*g2 + beta22(M,N);
dg3_dt = -epsilon*g3 + beta32(M,N);
dg4_dt = -epsilon*g4 + beta42(M,N);
system = [dg1_dt; dg2_dt; dg3_dt; dg4_dt];
In general there are no modifications needed after this line.
Differentiate and generate code (directional derivatives)¶
Exporting it to <system_name>.m
. This method uses directional derivatives.
Then using polarization identities derivatives can be calculated in arbitrary
direction.
suc = generate_directional_derivatives(...
system,... % n x 1 array of derivative symbolic expressions
coords,... % 1 x n array of symbols for states
par,... % 1 x np array of symbols used for parameters
system_name,... % argument specifying the system name
[matcontpath, 'Systems/']... % directory to save to file
);
Higher-order parameter-dependent multi-linear form.¶
Exporting it to <system_name>_multilinearforms.m
. These multi-linear forms are
currently only used in the computation of the parameter-dependent center
manifold for the codimension two Bogdanov-Takens bifurcation.
order = 3;
suc = generate_multilinear_forms(system_name, system, coords, par, order, ...
[matcontpath, 'Systems/']);