Fixed fundamental error in code for variance calculation

matlab_development/uq/examples/linear_harmonic_oscillator/demo.m

@@ -89,11 +89,12 @@ if not(coco_exist(out_run_name, 'run')) || true
     prob = coco_add_func(prob, 'orig.po', @linode_het_bc_v0, data, 'zero', ...
        'uidx', uidx([maps.x0_idx; maps.x1_idx; maps.T0_idx; maps.T_idx]));
 
-    prob = coco_add_pars(prob, 'xmax', maps.x0_idx(1), 'xmax', 'inactive');
+    prob = coco_add_pars(prob, 'xmax', maps.x0_idx, {'xmax', 'vmax'}, ...
+                         'inactive');
 
     % Vary the stiffness, k.
     bd = coco(prob, out_run_name, [], 1, ...
-              {'k', 'phi', 'xmax', 'om'}, [3,4]);
+              {'k', 'phi', 'xmax', 'vmax', 'om'}, [3,4]);
 else
     bd = coco_bd_read(out_run_name);
 end

matlab_development/uq/examples/linear_harmonic_oscillator/demo_uq.m

@@ -11,13 +11,67 @@ addpath(funcpath)
 clear all;
 %%
 run_name = 'init';
-if not(coco_exist(run_name, 'run'))
-    fprintf('Increase k to 1.0\n')
+if not(coco_exist(run_name, 'run')) 
+    fprintf('Increase k to 3.0\n')
     bd = demo(run_name);
 else
     bd = coco_bd_read(run_name);
 end
 
+last_run_name = run_name;
+run_name = 'sweep';
+if not(coco_exist(run_name, 'run'))
+    fprintf('Frequency Sweep for k=3\n') 
+    % Collect the problem instance from the last run
+    NTST = 60;
+    NCOL = 4;
+        
+
+    labs = coco_bd_labs(bd, 'EP');
+    lab_number = 1; % first entry
+    if lab_number < 0
+        num_labs = size(labs,2);
+        lab_number = num_labs + lab_number + 1;
+        lab = labs(lab_number);
+    else
+        lab = labs(lab_number);
+    end
+
+    sol = coll_read_solution('orig', last_run_name, lab);
+    % Reset Time
+    t0 = sol.tbp - sol.tbp(1);
+    x0 = sol.xbp;
+    p0 = sol.p;
+
+    prob = coco_prob(); % Reset run parameters
+    prob = coco_set(prob, 'ode', 'autonomous', false);
+    prob = coco_set(prob, 'coll', 'NTST', NTST, 'NCOL', NCOL);%, 'var', true);
+    prob = coco_set(prob, 'cont', 'PtMX', 1000, 'h', 0.5, 'almax', 30);    
+
+    coll_args = {@linode_het, @linode_het_DFDX, @linode_het_DFDP, ...
+      @linode_het_DFDT, t0, x0, {'k', 'phi', 'om'}, p0};
+
+    prob = ode_isol2coll(prob, 'orig', coll_args{:});
+
+    % Reapply Boundary Conditions
+    [data, uidx] = coco_get_func_data(prob, 'orig.coll', 'data', 'uidx');
+    maps = data.coll_seg.maps;
+
+    % _v0 boundary condition function locks in the initial velocity of zero
+    % which was found in the previous continuation run.
+    prob = coco_add_func(prob, 'orig.po', @linode_het_bc_v0, data, 'zero', ...
+       'uidx', uidx([maps.x0_idx; maps.x1_idx; maps.T0_idx; maps.T_idx]));
+
+    prob = coco_add_pars(prob, 'xmax', maps.x0_idx, {'xmax', 'vmax'}, ...
+                         'inactive');
+
+    % Vary the stiffness, k.
+    bd_sweep = coco(prob, run_name, [], 1, ...
+              {'om', 'phi', 'xmax', 'vmax', 'k'}, [0.01 10]);    
+else
+    bd_sweep = coco_bd_read(run_name);
+end
+
 %% Step 2
 %  Generate a sample points for the stochastic parameters.  In this code
 %  this will be done using the tensor product quadrature method of 
@@ -49,7 +103,7 @@ k_vals = mu_k + sig_k*nds;
 %  sample point.  Introduce zero functions that tie problem parameters
 %  together between problem instances.
 
-last_run_name = run_name;
+last_run_name = 'init';
 run_name = 'k_vals';
 if not(coco_exist(run_name, 'run'))
     % Collect the problem instance from the last run
@@ -287,7 +341,7 @@ end
 %     ps{i} = strcat('phi',int2str(i));
 % end
 
-bd = coco(prob, run_name, [], 1, [ps, {'mean', 'variance'}], [0.1, 10]);
+bd = coco(prob, run_name, [], 1, [ps, {'mean', 'variance'}], [0.01, 10]);
 %% 
 
 % rmpath(funcpath)
\ No newline at end of file

matlab_development/uq/toolbox/uq_pce_coefficients.m

@@ -2,6 +2,12 @@ function [data, y] = uq_pce_coefficients(prob, data, u)
 %UQ_PCE_COEFFICIENTS Summary of this function goes here
 %   Evaluate Response function and calculate coefficients of polynomial
 %   chaos expansion.
+%   u contains all Basis point and parameter values for M trajectories as
+%   well as continuation variables for the alphas where M is the numerical
+%   integration order of the PCE.
+%   u = [ xbp1, xbp2, ... xbpM, p1, p2, ... pM ,alpha];
+%
+%   y = alpha - Psi'*r(x,p);
 
 for i=1:data.M   % M is the number of quadrature points
      % Separate out components
@@ -20,11 +26,6 @@ for i=1:data.M   % M is the number of quadrature points
          r(i,:) = data.rhan(x, p);
      end
 end
-% last_state = data.M^data.s;
-% x = u(1:last_state);
-% alphas = u(last_state+1:end);
-% p = [];
-% r = data.fhan(x, p);
 
 alphas = u(data.pidx(end)+1:end);
 R = data.wtd_psi_mat'*r;

matlab_development/uq/toolbox/uq_pce_variance.m

@@ -5,7 +5,7 @@ function [data, y] = uq_pce_variance(prob, data, alphas)
 
 % For generic response functions that may have more than one output need to
 % provide data that allows appropriate reshaping of the array of alphas.
-y = sum(alphas(:, 1:end).^2);
+y = sum(alphas(:, 2:end).^2);
 
 end