Commiting toolbox and development files (instead of just .gitignore)

archive/demo.m

+% Parameter values
+% p(1) = Period, T
+% p(2) = Forcing Amplitude, A
+% p(3) = Linear Stiffness, k
+% p(4) = Damping, d
+% p(5) = Cubic Stiffness, eps
+p0 = [2*pi; 0.3; 1; 2*0.1; 5];
+u0 = [0; 1; 0];
+
+% Forward Simulation for Starting Guess for Periodic Orbit
+f = @(t,u) duff(u,p0);
+% Long Run Simulation to let transients die out
+[~, u0] = ode45(f, [0, 100*pi], u0);
+% Estimate for steady state solution
+[t0, u0] = ode45(f, [0, p0(1)], u0(end,:));
+
+%%
+% Initialize Boundary Conditions
+data = struct();
+data.fhan = @duff;
+data = per_bc_update(data, [], u0(1,:)', [], p0);
+
+funcs = { @duff @duff_DU @duff_DP };
+pnames = {'T' 'A' 'k' 'd' 'eps'};
+coll_args = {funcs{:} t0 u0 p0};
+bvp_args = { @per_bc, @per_Jbc, data, @per_bc_update };
+
+%%
+prob = coco_prob();
+prob = coco_set(prob, 'coll', 'NTST', 15, 'var', true);
+% prob = coco_set(prob, 'cont', 'NAdapt', 1);
+prob = ode_isol2bvp(prob, '', coll_args{:}, pnames, bvp_args{:});
+
+%%
+[data, uidx] = coco_get_func_data(prob, 'bvp.seg1.coll', 'data', 'uidx');
+
+maps = data.coll_seg.maps;
+
+prob = coco_add_glue(prob, 'glue', uidx(maps.T_idx), uidx(maps.p_idx(1)));
+prob = coco_add_pars(prob, 'x_max', uidx(maps.x0_idx(1)), {'x_max'});
+
+%%
+bd = coco(prob, 'optim', [], 1, {'T','bvp.seg1.coll.T0','x_max'}, [2*pi/1.3, 2*pi/0.7]);
+

archive/demo_coll.m

+% Parameter values
+% p(1) = frequency of oscillation, w
+% p(2) = Forcing Amplitude, A
+% p(3) = Linear Stiffness, k
+% p(4) = Damping, d
+% p(5) = Cubic Stiffness, eps
+% p0 = [2*pi; 0.3; 1; 2*0.1; 5];
+p0 = [1; 0; 1; 0; 0];
+% u0 = [0; 1; 0; 1];
+
+t0 = linspace(0,1,200)';
+u0 = [cos(t0), -sin(t0), sin(t0), cos(t0)];
+
+% % Forward Simulation for Starting Guess for Periodic Orbit
+% f = @(t,u) duff(u,p0);
+% % Long Run Simulation to let transients die out
+% [t0, u0] = ode45(f, [0, 2*pi], u0);
+
+% Estimate for steady state solution
+% [t0, u0] = ode45(f, [0, 2*pi/p0(1)], u0(end,:));
+data = struct();
+data.fhan = @duff;
+data = per_bc_update(data, [], u0(1,:)', [], p0);
+
+prob = coco_prob();
+funcs = { @duff @duff_DU @duff_DP };
+pnames = {'w' 'A' 'k' 'd' 'eps'};
+
+coll_args = {funcs{:}, t0, u0, p0};
+
+bvp_args  = {@per_bc, @per_Jbc, data, @per_bc_update};
+prob = ode_isol2bvp(prob, '', coll_args{:}, pnames, bvp_args{:});
+
+[data, uidx] = coco_get_func_data(prob, 'bvp.seg1.coll', 'data', 'uidx');
+maps = data.coll_seg.maps;
+
+prob = coco_add_glue(prob, 'glue', uidx(maps.T_idx), uidx(maps.p_idx(1)));
+prob = coco_add_glue(prob, 'glue2', uidx(maps.T0_idx), uidx(maps.x0_idx(3)));
+prob = coco_add_pars(prob, 'v0', uidx(maps.x0_idx(2)), 'v0', 'active');
+prob = coco_add_pars(prob, 'x_max', uidx(maps.x0_idx(1)), 'x_max', 'active');
+
+% Start continuation of periodic orbits with updates to the Poincare
+% section before each continuation step.
+
+prob = coco_set(prob, 'cont', 'NAdapt', 10);
+cont_args = { 0, {'w', }};
+
+coco(prob, 'optim', [], cont_args{:});
\ No newline at end of file

archive/demo_po.m

+% Parameter values
+% p(1) = Period, T
+% p(2) = Forcing Amplitude, A
+% p(3) = Linear Stiffness, k
+% p(4) = Damping, d
+% p(5) = Cubic Stiffness, eps
+p0 = [2*pi; 0.3; 1; 2*0.1; 5];
+u0 = [0; 1];
+
+% Forward Simulation for Starting Guess for Periodic Orbit
+f = @(t,u) duffna(t,u,p0);
+
+% Long Run Simulation to let transients die out
+[~, u0] = ode45(f, [0, 100*pi], u0);
+
+% Estimate for steady state solution
+[t0, u0] = ode45(f, [0, p0(1)], u0(end,:)');
+
+funcs = { @duffna @duffna_DU @duffna_DP @duffna_DT};
+pnames = {'T' 'A' 'k' 'd' 'eps'};
+coll_args = {funcs{:} t0 u0 pnames p0};
+prob = coco_prob();
+prob = coco_set(prob, 'ode', 'autonomous', false);
+prob = coco_set(prob, 'coll', 'NTST', 15);
+prob = ode_isol2po(prob, '', coll_args{:});
+
+% Adaptive discretization before each continuation step
+prob = coco_set(prob, 'cont', 'NAdapt', 1);
+
+[data, uidx] = coco_get_func_data(prob, 'po.orb.coll', 'data', 'uidx');
+maps = data.coll_seg.maps;
+
+prob = coco_add_glue(prob, 'glue', uidx(maps.T_idx), uidx(maps.p_idx(1)));
+prob = coco_add_pars(prob, 'v0', uidx(maps.x0_idx(2)), 'v0', 'active');
+
+
+prob = coco_add_pars(prob, 'x_max', uidx(maps.x0_idx(1)), {'x_max'});
+
+bd = coco(prob, 'optim', [], 1, {'po.period', 'T', 'x_max'}, [pi, 4*pi]);
+bd2 = coco_bd_read('optim');
+u = coco_bd_col(bd2, {'T', '||po.orb.x||'});
+plot(2*pi./u(1,:),u(2,:))
+

archive/demo_po_a.m

+% Parameter values
+% p(1) = Period, T
+% p(2) = Forcing Amplitude, A
+% p(3) = Linear Stiffness, k
+% p(4) = Damping, d
+% p(5) = Cubic Stiffness, eps
+p0 = [2*pi; 0.3; 1; 2*0.1; 5];
+u0 = [0; 1; 0];
+
+% Forward Simulation for Starting Guess for Periodic Orbit
+f = @(t,u) duff(u,p0);
+
+% Long Run Simulation to let transients die out
+[~, u0] = ode45(f, [0, 20*pi], u0);
+
+% Estimate for steady state solution
+[t0, u0] = ode45(f, [0, p0(1)], u0(end,:)');
+
+funcs = {@duff @duff_DU @duff_DP};
+pnames = {'T' 'A' 'k' 'd' 'eps'};
+coll_args = {funcs{:} t0 u0 pnames p0};
+
+prob = ode_isol2po(coco_prob(), '', coll_args{:});
+prob = coco_set(prob, 'coll', 'NTST', 15);
+
+% Adaptive discretization before each continuation step
+prob = coco_set(prob, 'cont', 'NAdapt', 1);
+
+[data, uidx] = coco_get_func_data(prob, 'po.orb.coll', 'data', 'uidx');
+maps = data.coll_seg.maps;
+
+prob = coco_add_glue(prob, 'glue', uidx(maps.T_idx), uidx(maps.p_idx(1)));
+prob = coco_add_pars(prob, 'v0', uidx(maps.x0_idx(2)), 'v0', 'active');
+
+
+prob = coco_add_pars(prob, 'x_max', uidx(maps.x0_idx(1)), {'x_max'});
+
+bd = coco(prob, 'optim', [], 1, {'po.period', 'T', 'x_max'}, [pi, 4*pi]);
+bd2 = coco_bd_read('optim');
+u = coco_bd_col(bd2, {'T', '||po.orb.x||'});
+plot(2*pi./u(1,:),u(2,:))
+

archive/duff.m

+function y = duff(u,p)
+
+% Autonomous, Vectorized encoding of the Duffing oscillator
+% x'' + d*x' + k*x + eps*x^3 = A*cos(2*pi*t/T)
+% u1 = x
+% u2 = x'
+% u3 = sin(wt)
+% u4 = cos(wt)
+
+w = p(1,:);
+A = p(2,:);
+k = p(3,:);
+d = p(4,:);
+eps = p(5,:);
+
+u1 = u(1,:); % Position
+u2 = u(2,:); % Velocity
+u3 = u(3,:); % sin(wt)
+u4 = u(4,:); % cos(wt)
+
+y(1,:) = u2;
+y(2,:) = A.*u4 - d.*u2- k.*u1 - eps.*u1.^3;
+% Solution to u3 and u4 is u3 = sin(wt), u4 = cos(wt)
+% Therefore, this is an autonomous way to encode the non-autonomous forcing
+% function.  See Mehdi RSPA paper, reference ~47 (the one about
+% continuation in a finite element model)
+y(3,:) = u3 + w.*u4 - u3.*(u3.^2 + u4.^2);
+y(4,:) = u4 - w.*u3 - u4.*(u3.^2 + u4.^2);
+
+end
\ No newline at end of file

archive/duff_DP.m

+function J = duff_DP(u,p)
+
+u1 = u(1,:);
+u2 = u(2,:);
+u3 = u(3,:);
+u4 = u(4,:);
+
+J = zeros(4,size(p,1),size(u,2));
+
+J(2,2,:) = u4;
+J(2,3,:) = -u1;
+J(2,4,:) = -u2;
+J(2,5,:) = -u1.^3;
+J(3,1,:) = u4;
+J(4,1,:) = -u3;
+
+end
\ No newline at end of file

archive/duff_DU.m

+function J = duff_DU(u,p)
+
+u1 = u(1,:);
+u3 = u(3,:);
+u4 = u(4,:);
+w = p(1,:);
+A = p(2,:);
+k = p(3,:);
+d = p(4,:);
+eps = p(5,:);
+
+J = zeros(4,4,numel(u1));
+
+J(1,2,:) = 1;
+J(2,1,:) = -k - 3*eps.*u1.^2;
+J(2,2,:) = -d;
+J(2,4,:) = A;
+J(3,3,:) = 1 - 3*u3.^2 - u4.^2;
+J(3,4,:) = w - 2*u3.*u4;
+J(4,3,:) = -w - 2*u3.*u4;
+J(4,4,:) = 1 - u3.^2 - 3*u4.^2;
+
+end
\ No newline at end of file

archive/duffna.m

+function y = duffna(t,u,p)
+
+% non-Autonomous, Vectorized encoding of the Duffing oscillator
+% x'' + d*x' + k*x + eps*x^3 = A*cos(2*pi*t/T)
+% u1 = x
+% u2 = x'
+
+T = p(1,:);
+A = p(2,:);
+k = p(3,:);
+d = p(4,:);
+eps = p(5,:);
+
+u1 = u(1,:); % Position
+u2 = u(2,:); % Velocity
+
+y(1,:) = u2;
+y(2,:) = A.*cos((2*pi./T).*t) - d.*u2- k.*u1 - eps.*u1.^3;
+
+end
\ No newline at end of file

archive/duffna_DP.m

+function J = duffna_DP(t,u,p)
+
+T = p(1,:);
+A = p(2,:);
+
+u1 = u(1,:);
+u2 = u(2,:);
+
+J = zeros(2,numel(p(:,1)),numel(u1));
+
+J(2,1,:) = (2*pi*A.*t./(T.^2)).*sin((2*pi./T).*t);
+J(2,2,:) = cos(2*pi./T.*t);
+J(2,3,:) = -u1;
+J(2,4,:) = -u2;
+J(2,5,:) = -u1.^3;
+
+end
\ No newline at end of file

archive/duffna_DT.m

+function Jt = duffna_DT(t,u,p)
+
+T = p(1,:);
+A = p(2,:);
+
+Jt = zeros(2,numel(t));
+Jt(2,:) = -(2*pi./T).*A.*sin((2*pi./T).*t);
+
+end
\ No newline at end of file

archive/duffna_DU.m

+function J = duffna_DU(t,u,p)
+
+u1 = u(1,:);
+T = p(1,:);
+A = p(2,:);
+k = p(3,:);
+d = p(4,:);
+eps = p(5,:);
+
+J = zeros(2,2,numel(u1));
+
+J(1,2,:) = 1;
+J(2,1,:) = -k - 3*eps.*u1.^2;
+J(2,2,:) = -d;
+
+end
\ No newline at end of file

archive/per_Jbc.m

+function  Jbc = per_Jbc(data, T, x0, x1, p)
+
+Jbc = data.J;
+
+
+end
\ No newline at end of file

archive/per_bc.m

+function [data, y] = per_bc(T, u0, u1, p)
+
+% Boundary Conditions for the oscillator
+% 1. u1(0) - u1(T) = 0 -----| Initial Postion = Final Position
+% 2. u2(0) - u2(T) = 0 -----| Initial Velocity = Final Velocity
+% 3. u3(0) - u3(1) = 0 -----| 
+% 4. u4(0) - u4(1) = 0 -----|
+% 4. u2(0) = 0 -------------| Phase Condition: Start orbit where the Velocity 
+%                           |                  is zero (max Amplitude).
+
+u01 = u0(1,:); % Initial Position
+u02 = u0(2,:); % Initial Velocity
+u03 = u0(3,:);
+u04 = u0(4,:);
+
+u11 = u1(1,:); % Final Position
+u12 = u1(2,:); % Final Velocity
+u13 = u1(3,:); 
+u14 = u1(4,:); 
+
+y(1,:) = u01 - u11;       % Periodic Position BC
+y(2,:) = u02 - u12;       % Periodic Velocity BC
+y(3,:) = u03 - u13;
+y(4,:) = u04 - u14;
+y(5,:) = u02;             % Poincare Section
+
+end
\ No newline at end of file

archive/per_bc_update.m

+function data = per_bc_update(data, T, u0, u1, p)
+
+% Update Function for the boundary conditions
+
+n = numel(u0);
+q = numel(p);
+
+pmat = sparse(n,q);
+pmat(n,1) = 1;
+
+% Jacobian for Poincare Section at u2 = 0
+data.J = [ sparse(n,1), speye(n,n),                -speye(n,n),  pmat;
+           sparse(1,1), sparse([0 1 zeros(1,n-2)]), sparse(1,n), sparse(1,q)];
+
+end
\ No newline at end of file

coco_example_runs/duff/F.m

+function y = F(u,p)
+
+% State Vector Assignment
+u1 = u(1,:);        % Position
+u2 = u(2,:);        % Velocity
+u3 = u(3,:);        % time
+
+% Parameter Assignment
+A = p(1,:);         % Forcing function amplitude
+k = p(2,:);         % Linear Stiffness
+d = p(3,:);         % Linear Damping
+w = p(4,:);         % Forcing function frequency
+eps = p(5,:);       % Nonlinear Stiffness
+
+y(1,:) = u2;
+y(2,:) = A.*cos(w.*u3) - d.*u2-k.*u1 - eps.*u1.^3;
+y(3,:) = 1;
+
+end
\ No newline at end of file

coco_example_runs/duff/F_dP.m

+function J = F_dP(u,p)
+
+% State Vector Assignment
+u1 = u(1,:);        % Position
+u2 = u(2,:);        % Velocity
+u3 = u(3,:);        % time
+
+A = p(1,:);
+w = p(4,:);
+
+J = zeros(3,5,numel(u1));
+J(2,1,:) = cos(w.*u3);
+J(2,2,:) = -u1;
+J(2,3,:) = -u2;
+J(2,4,:) = -A.*u3.*sin(w.*u3);
+J(2,5,:) = -u1.^3;
+
+end
\ No newline at end of file

coco_example_runs/duff/F_dU.m

+function J = F_dU(u,p)
+
+% State Vector Assignment
+u1 = u(1,:);        % Position
+u3 = u(3,:);        % Position
+
+% Parameter Assignment
+A = p(1,:);         % Forcing function amplitude
+k = p(2,:);         % Linear Stiffness
+d = p(3,:);         % Linear Damping
+w = p(4,:);         % Forcing function frequency
+eps = p(5,:);       % Nonlinear Stiffness
+
+% Assign Jacobian Terms
+J = zeros(3,3,numel(u1));
+
+J(1,2,:) = 1;
+J(2,1,:) = -k - 3*eps.*u1.^2;
+J(2,2,:) = -d;
+J(2,3,:) = -A.*w.*sin(w.*u3);
+
+end
\ No newline at end of file

coco_example_runs/duff/adj_F.m

+function [data, y] = adj_F(prob, data, u)
+
+[data, uidx] = coco_get_func_data(prob, 'coll', 'data', 'uidx');
+
+
+end
\ No newline at end of file