updated comments

mains/confusion.m

@@ -0,0 +1,39 @@
+% cannot explain why following is stable (set rho to pi/8)
+
+init; [ts,ys] = EPBM_PMFCmS(LF(pars), NF, 2 * tspan, y0, 2 * 800, struct('z', [-1 ; legpts(4)], 'b', -1 * ones(5), 'mS', 1, 'kappa', 2, 'parameters', pars)); mesh(sfilter(ys))
+init; [ts,ys] = EPBM_PMFCmS(LF(pars), NF, 2 * tspan, y0, 2 * 800, struct('z', lobpts(4), 'b', -1 * ones(4), 'mS', [], 'kappa', 1, 'parameters', pars)); mesh(sfilter(ys))
+
+% using rho = pi/64, and tspan = [0 80]
+init; [ts,ys] = EPBM_PMFCmS(LF(pars), NF, 2 * tspan, y0, 4 * 800, struct('z', [-1 ; legpts(4)], 'b', -1 * ones(5), 'mS', 1, 'kappa', 0, 'parameters', pars)); mesh(sfilter(ys))
+% however, it gets worse if kappa > 0?
+
+
+
+%%
+
+
+z1s = linspace(0, 20,200);
+z2s = linspace(-2,2,100);
+a = zeros(length(z1s), length(z2s));
+a_new = zeros(length(z1s), length(z2s));
+rho0 = pi / 8;
+
+for i = 1 : length(z1s)
+    shift = z1s(i)^(3/3) / tan(pi/2 + rho0);
+    a(i,:) = ampETDPBM(1i * z1s(i) + shift, 1i * z2s(:) - shift, struct('z', lobpts(4), 'b', -1 * ones(4), 'mS', [], 'kappa', 3, 'alpha', 2));%struct('z', [-1 ; legpts(4)], 'b', -1 * ones(5), 'mS', 1, 'kappa', 2, 'alpha', 2)); 
+    %a(i,:) = ampETDPBM(1i * z1s(i) + shift, 1i * z2s(:) - shift, struct('z', [-1 ; legpts(4)], 'b', -1 * ones(5), 'mS', 1, 'kappa', 1, 'alpha', 1)); 
+%    a(i,:) = abs(rERK4(1i * z1s(i) + shift, 1i * z2s(j) - shift));
+%     for j = 1 :length(z2s)
+%         shift = z1s(i) / tan(pi/2 + rho0);
+%         a(i,j) = ampETDPBM(1i * z1s(i) + shift, 1i * z2s(j) - shift, struct('z', [-1 ; legpts(4)], 'b', -1 * ones(5), 'mS', 1, 'kappa', 0, 'alpha', 2));
+%         %a(i,j) = abs(rERK4(1i * z1s(i) + shift, 1i * z2s(j) - shift));
+%     end
+end
+
+figure(100);
+
+%a(a > 1) = NaN;
+surf(z2s, z1s, a); shading interp; view([90 90]); hold on;
+[~, cnt] = contour(z2s, z1s, a, [0 1], 'k', 'LineWidth', 2); hold off;
+cnt.ZLocation = 1;
+

mains/experiment/trash/main_epsilon_graph.m

@@ -0,0 +1,213 @@
+% include directories
+addpath(genpath('../../stepper/common'));
+addpath(genpath('../../stepper/rk'));
+addpath(genpath('../../stepper/poly'));
+addpath(genpath('../../stepper/sdc'));
+
+%% == load equation =======================================================
+equation = 'zds';
+
+wd = pwd();
+cd(fullfile('../../stepper/equations/', equation))
+init
+cd(wd);
+
+%% experiment settings
+integrator = 'erk'; % erk, esdc, epbm
+partitioning = 'hyperdiffusion-8'; % rotation, rotation-uxx, translation, hyperdiffusion
+
+par_to_label_value = @(v) num2str(v);
+par_to_label_name  = @(v) ['\', v];
+switch( partitioning )
+    case 'rotation'
+        par_name = 'rho';
+        par_values = [0 pi/256 pi/128 pi/64 pi/32];  
+        par_to_label_value = @rhoLabel; 
+    case 'rotation-uxx'
+        par_name = 'rhouxx';
+        par_values = [0 pi/256 pi/128 pi/64 pi/32];
+        par_to_label_value = @rhoLabel;
+    case 'translation'
+        par_name = 'epsilon';
+        par_values = [0 1 2 4 8];
+        par_to_label_value = @epsilonLabel;
+    case {'hyperdiffusion-2', 'hyperdiffusion-4', 'hyperdiffusion-6', 'hyperdiffusion-8'}
+
+        par_name   = 'hyperv_coeff';
+        par_to_label_name  = @(v) '\omega';
+        par_to_label_value = @omegaLabel;
+
+        switch(partitioning)
+            case 'hyperdiffusion-2'
+                 pars.hyperv_order = 2;
+                par_values = [0 1e6 1e8 1e10 1e12];
+            case 'hyperdiffusion-4'
+                pars.hyperv_order = 4;
+                par_values = [0 1e4 1e6 1e8 1e10];
+            case 'hyperdiffusion-6'
+                pars.hyperv_order = 6;
+                par_values = [0 1e2 1e4 1e6 1e8];
+            case 'hyperdiffusion-8'
+                pars.hyperv_order = 8;
+                par_values = [0 1 1e2 1e4 1e6];
+        end
+
+    otherwise
+        error('invalid partitioning');
+end
+
+switch ( integrator )
+    case 'erk'
+        integrator = @ERK;
+        integrator_options = struct('coeffGenerator', @ERK4, 'max_ts_to_store', 2);
+        line_marker = 's-';
+    case 'esdc'
+        integrator = @ESDC;
+        integrator_options = struct('tau', lobpts(4, [0 1]), 'm', 6); 
+        line_marker = 'd-';
+    case 'epbm'
+        integrator = @EPBM_PMFCmS;
+        integrator_options = struct('z', [-1 ; legpts(4)], 'b', -1 * ones(5), 'mS', 1, 'kappa', 0, 'alpha', 1);
+        %integrator_options = struct('z', [-1 ; legpts(4)], 'b', -1 * ones(5), 'mS', 1, 'kappa', 1, 'alpha', 1);
+        %integrator_options = struct('z', [-1 ; legpts(4)], 'b', -1 * ones(5), 'mS', 1, 'kappa', 2, 'alpha', 2);
+        %integrator_options = struct('z', lobpts(4), 'b', -1 * ones(4), 'mS', [], 'kappa', 1, 'alpha', 2);
+        line_marker = '*-';
+    otherwise
+        error('invalid integrator');
+end
+
+%% == compute reference using IMEXRK4 =====================================
+fprintf('Computing reference... ');
+reference_ppars = mergeStructs(struct('delta', 0, 'epsilon', 0), pars, false);
+reference_options = struct('parameters', reference_ppars);
+L_ref = LF(reference_ppars);
+RHS = @(t,y,pars) L_ref .* y + NF(t,y,pars);
+[~, y_reference] = rk4(RHS, tspan, y0(:), Nt_reference, reference_options);
+% reference_options = struct('coeffGenerator', @IMRK4, 'parameters', reference_ppars);
+% [~, ys] = IMRK(LF(reference_ppars), NF, tspan, y0, Nt_reference, reference_options);
+% y_reference = ys(:,end);
+fprintf('done.\n');
+
+%% == prepare run parameters ==============================================
+num_pars = length(par_values);
+runs     = cell(num_pars, 1);
+for i = 1 : num_pars
+    prb_pars = mergeStructs(struct(par_name, par_values(i)), pars, false); % problem parameters
+    runs{i} = struct( ...
+        'integrator', integrator, ...
+        'options', mergeStructs(struct('parameters', prb_pars), integrator_options), ...
+        'legend', ['$', par_to_label_name(par_name), ' = ', par_to_label_value(par_values(i)), '$']);
+end
+
+% -- add IMRK4 integrator ---
+imrk4_options = struct('coeffGenerator', @IMRK4, 'parameters', reference_ppars);
+runs{end+1} = struct('integrator', @IMRK, 'options', imrk4_options, 'legend', 'IMRK4');
+
+%% == run integrators =====================================================
+num_Nts  = length(Nts);
+num_runs = length(runs);
+errors = zeros(num_Nts, num_runs);
+parfor i = 1 : num_runs
+    prob_pars = runs{i}.options.parameters
+    fprintf(' Running %s = %i\n',  par_name, prob_pars.(par_name));
+    for j = 1 : num_Nts
+
+        prob_pars.method_stepsize = tspan(end) / Nts(j);
+        fprintf('     Nt = %i/%i...', j , num_Nts);
+
+        [ts, ys] = runs{i}.integrator(LF(prob_pars), NF, tspan, y0, Nts(j), runs{i}.options);
+        errors(j, i) = error_norm(ys(:,end), y_reference);
+
+        fprintf('done.\n');
+    end
+end
+
+%% == produce plots (OLD) =======================================================
+% f = figure();
+% 
+% hs = diff(tspan) ./ Nts;
+% legend_labels = cellfun(@(run) run.legend, runs, 'UniformOutput', false);
+% 
+% h = loglog(hs, errors(:, 1:end-1), line_marker, 'LineWidth', 3, 'MarkerSize', 10);
+% set(h, {'MarkerFaceColor'}, get(h, 'Color'))
+% xlabel('Stepsize $h$', 'Interpreter', 'Latex');
+% ylabel('Relative Error','Interpreter', 'Latex');
+% axis([hs(end) hs(1) 1e-7 1e2])
+% 
+% hold on;
+% loglog(hs, errors(:, end), 'k--', 'LineWidth', 2);
+% hold off;
+% 
+% legend(legend_labels, 'Location', 'NorthWest', 'Interpreter', 'latex'); legend boxoff;
+% 
+% f_info = functions(integrator);
+% exportFigure(f, struct('Format', 'pdf', 'SavePath', fullfile('output', [equation, '-', num2str(Nx), '-', f_info.function, '-', partitioning])));
+
+%% == produce plots =======================================================
+f = figure();
+
+hs = diff(tspan) ./ Nts;
+line_markers = {'-', '-<', '->', '-^', '-v'};
+
+MS = 10 - 5 * (0:num_pars - 1) / (num_pars - 1); MS = [MS(end) MS(1:end-1)];
+LW = 3.5 - 1 * (0:num_pars - 1) / (num_pars - 1); LW = [LW(end) LW(1:end-1)];
+
+for i = 1 : length(par_values)
+    h = loglog(hs, errors(:, i), line_markers{i}, 'LineWidth', LW(i), 'MarkerSize', MS(i)); hold on;
+    %h = loglog(hs, errors(:, i), line_marker, 'LineWidth', LW(i), 'MarkerSize', 8.75, 'MarkerIndices', i:num_pars:length(errors(:, i))); hold on;
+    set(h, 'MarkerFaceColor', min(1, get(h, 'Color') * 1))
+    set(h, 'MarkerEdgeColor', min(1, get(h, 'Color') * 1))
+end
+hold off;
+
+
+xlabel('Stepsize $h$', 'Interpreter', 'Latex');
+ylabel('Relative Error','Interpreter', 'Latex');
+axis([hs(end) hs(1) 1e-8 1e2])
+
+% place rho = 0 line on top.
+chH = get(gca,'Children');
+set(gca,'Children',[chH(end); chH(1:end-1)]);
+
+hold on;
+loglog(hs, errors(:, end), 'k--', 'LineWidth', 2);
+hold off;
+
+xticks([1e-3, 1e-2, 1e-1])
+yticks(10.^(-8:2:2))
+%legend(legend_labels, 'Location', 'NorthWest', 'Interpreter', 'latex'); legend boxoff;
+
+f_info = functions(integrator);
+exportFigure(f, struct('Format', 'pdf', 'SavePath', fullfile('output', [equation, '-', num2str(Nx), '-', f_info.function, '-', partitioning]), 'PaperPosition', [0 0 8 7]));
+
+
+%% Legend Figure
+f = figure();
+
+legend_labels = cellfun(@(run) run.legend, runs, 'UniformOutput', false);
+for i = 1 : length(par_values)
+    h = plot(0, 0, line_markers{i}, 'LineWidth', LW(i), 'MarkerSize', MS(i));
+    set(h, 'MarkerFaceColor', min(1, get(h, 'Color') * 1))
+    set(h, 'MarkerEdgeColor', min(1, get(h, 'Color') * 1))
+    hold on;
+end
+legend(legend_labels, 'Location', 'northoutside', 'Interpreter', 'latex', 'Orientation', 'horizontal'); legend boxoff;
+set(gca, 'FontSize', 12);
+exportFigure(f, struct('Format', 'pdf', 'SavePath', fullfile('output', ['legend-', partitioning]), 'PaperPosition', [0 0 20 8]));
+
+function l = rhoLabel(rho)
+    if(rho == 0)
+        l = '0 \hspace{1em}';
+    else
+        [~, b] = rat(rho / pi);
+        l = ['\pi / ', num2str(b), '\hspace{1em}'];
+    end
+end
+
+function l = epsilonLabel(epsilon)
+    l = [num2str(epsilon), '\hspace{1em}'];
+end
+
+function l = omegaLabel(omega)
+    l = [num2str(omega), '\hspace{1em}'];
+end

mains/experiment/trash/main_epsilon_graph_old.m

@@ -0,0 +1,89 @@
+% include directories
+addpath(genpath('../../stepper/common'));
+addpath(genpath('../../stepper/rk'));
+
+% experiment settings
+equation = 'zds';
+epsilons = [0 0.5 1 2 4 6 8 10];
+deltas = linspace(0, 8/100, 9);
+integrator = @ERK;
+integrator_params = struct('coeffGenerator', @ERK4, 'max_ts_to_store', 2);
+
+%% == load equation =======================================================
+wd = pwd();
+cd(fullfile('../../stepper/equations/', equation))
+init
+cd(wd);
+
+%% == compute reference using IMEXRK4 =====================================
+fprintf('Computing reference... ');
+pars.epsilon = 0;
+[~, ys] = IMRK(LF(pars), NF, tspan, y0, Nt_reference, struct('coeffGenerator', @IMRK4, 'parameters', pars));
+y_reference = ys(:,end);
+fprintf('done.\n');
+
+%% == run experiment with reference integrator ============================
+num_Nts   = length(Nts);
+imrk_error = zeros(num_Nts, 1);
+
+pars.epsilon = 0; pars.delta = 0;
+fprintf(' Running Reference Integrator = %i\n',  pars.epsilon);
+for j = 1 : num_Nts
+    fprintf('     Nt = %i/%i...', j , num_Nts);
+    [~, ys] = IMRK(LF(pars), NF, tspan, y0, Nts(j), struct('coeffGenerator', @IMRK4, 'parameters', pars));
+    imrk_error(j) = error_norm(ys(:,end), y_reference);
+    fprintf('done.\n');
+end
+
+%% == run experiment with exponential integrator ==========================
+% num_epsilons = length(epsilons);
+% num_Nts = length(Nts);
+% errors = zeros(num_epsilons, num_Nts);
+% legend_labels = cell(num_epsilons, 1);
+% 
+% parfor i = 1 : num_epsilons
+%     epars = mergeStructs(struct('epsilon', epsilons(i)), pars, false);
+%     fprintf(' Running Epsilon = %i\n',  epars.epsilon);
+%     for j = 1 : num_Nts
+%         fprintf('     Nt = %i/%i...', j , num_Nts);
+%         
+%         [ts, ys] = integrator(LF(epars), NF, tspan, y0, Nts(j), mergeStructs(struct('parameters', epars), integrator_params, false));
+%         errors(i,j) = error_norm(ys(:,end), y_reference);
+%         
+%         fprintf('done.\n');
+%     end
+%     legend_labels{i} = ['\epsilon = ', num2str(epsilons(i))];
+% end
+
+num_deltas = length(deltas);
+num_Nts = length(Nts);
+errors = zeros(num_deltas, num_Nts);
+legend_labels = cell(num_deltas, 1);
+
+parfor i = 1 : num_deltas
+    epars = mergeStructs(struct('delta', deltas(i)), pars, false);
+    fprintf(' Running delta = %i\n',  epars.epsilon);
+    for j = 1 : num_Nts
+        fprintf('     Nt = %i/%i...', j , num_Nts);
+
+        [ts, ys] = integrator(LF(epars), NF, tspan, y0, Nts(j), mergeStructs(struct('parameters', epars), integrator_params, false));
+        errors(i,j) = error_norm(ys(:,end), y_reference);
+
+        fprintf('done.\n');
+    end
+    legend_labels{i} = ['\delta = ', num2str(deltas(i))];
+end
+
+%% == plot ================================================================
+hs = diff(tspan) ./ Nts;
+
+loglog(hs, errors);
+legend(legend_labels);
+xlabel('h');
+ylabel('error');
+axis([hs(end) hs(1) 1e-16 1e2])
+
+hold on;
+loglog(hs, imrk_error, 'k--');
+hold off;
+