Skip to content
Snippets Groups Projects
Commit 3954b48d authored by jcandrsn's avatar jcandrsn
Browse files

Updating .gitignore

parent 20e18c61
No related branches found
No related tags found
1 merge request!1Development
Showing
with 4 additions and 2184 deletions
......@@ -23,9 +23,10 @@ octave-workspace
# coco solution files
sol*.mat
coco_log.txt
coco_scr.txt
bd.mat
coco_log*.txt
coco_scr*.txt
bd*.mat
matlab_development/uq/examples/**/data
# Some Mac file
.DS_Store
......
% 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]);
% 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
% 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,:))
% 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,:))
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
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
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
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
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
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
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
This diff is collapsed.
function Jbc = per_Jbc(data, T, x0, x1, p)
Jbc = data.J;
end
\ No newline at end of file
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
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
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
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
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
File deleted
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment