Reworked predator prey filippov system

toolbox/examples/predatorpreyFilippov/generateEulerSol.m

@@ -0,0 +1,99 @@
+%% SETUP
+% parameter values p = (r1, r2, beta1, beta2, q1, q2, m, e)
+% order as in paper (see rhs file)
+m     = 0.790;
+r1    = 0.836;
+e     = 0.948;
+q1    = 0.772;
+aq    = 0.660;
+beta2 = 0.896;
+
+beta1 = 7.81; q2 = 1.5; r2 = 0.3;
+
+% initial values, parameters, timespan
+a = 0.286975;
+x0_1 = [a ; a ; r1-r2];
+p = [r1, r2, beta1, beta2, q1, q2, m, e, aq];
+tspan = [0 100];
+
+% configure plotting
+X_plot = linspace(tspan(1), tspan(end), 1000);
+fignum = 1002;
+figure(fignum); clf; hold('on');
+plotit = @plotter;
+
+% solver selection and configuration
+intEuler       = @explEuler;
+eulerStep      = 1e-7;
+namePlainEuler = @(f) sprintf('plain %s' , func2str(f));
+
+%% COMPUTATION
+% Now let the user decide if an euler solution is generated or loaded
+fprintf('\nThis is the Euler solution generation script. Proceed with generation?\n');
+choices = {"no (load from .mat file)", "yes (can take up to 20 minutes)"};
+default_choice_index = 1;
+[idx, val] = userchoice(choices, default_choice_index);
+doEuler = false;
+if idx == 2
+    doEuler = true;
+end
+
+% EULER Integration
+[owndir, ~] = fileparts(mfilename('fullpath'));
+euler_fname = fullfile(owndir, sprintf('sol_euler_red_%.0e.mat', eulerStep));
+EulerFileIsPresent = isfile(euler_fname);
+if EulerFileIsPresent && ~doEuler
+    fprintf('Loading sol_euler from file %s\n', euler_fname);
+    tmp = load(euler_fname, 'sol_euler_ds');
+    sol_euler = tmp.sol_euler_ds;
+    doEuler = true;
+else
+    if ~EulerFileIsPresent
+        disp("Euler solution file missing... Generating file");
+    else
+        disp("Computing Euler solution");
+    end
+    % Generate euler solution
+    fprintf('Integrating with integrator %s (might take a while) ...\n', func2str(intEuler))
+    figure(fignum);
+    th = tic();
+    sol_euler = intEuler(@(t,x) pprhs(t,x,p), tspan, x0_1, eulerStep);
+    time_euler = toc(th); fprintf('Euler took %g s\n', time_euler);
+    fprintf('Saving result to %s for later reuse.\n', euler_fname);
+    % reduce euler solution to make file smaller and save it
+    ds = 100;
+    idx = 1:ds:numel(sol_euler.x);
+    sol_euler_ds.x = sol_euler.x(idx);
+    sol_euler_ds.y = sol_euler.y(:, idx);
+    save(euler_fname, "sol_euler_ds");
+    doEuler = true; % in case we ended up here because the file did not exist
+end
+
+if doEuler
+    X_euler = X_plot;
+    Y_euler = transpose(interp1(sol_euler.x, transpose(sol_euler.y), X_euler));
+    linewidth = 2.0;
+    hEuler = plotit(fignum, Y_euler, 'c', namePlainEuler(intEuler), linewidth);
+end
+
+
+
+% FINITO
+return
+
+%% HELPERS
+
+function h = plotter(fignum, y, color, name, lw)
+   figure(fignum); hold on;
+   h = plot3(y(3,:), y(2,:), y(1,:), 'Color', color, 'LineWidth', lw, 'DisplayName', name);
+   view([97 51]);
+   grid on;
+   box on;
+   xlabel('Predator');
+   ylabel('Prey 2');
+   zlabel('Prey 1');
+   legend('location', 'northeast');
+   drawnow
+   pause(1.0);
+   set(fignum, 'Position', [200  250  750  375]);
+end

toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m → toolbox/examples/predatorpreyFilippov/predatorPrey3D_bench.m

@@ -44,6 +44,12 @@
 wrt_y = 3;
 eulerDisturbH = 1e-6;
 
+%% FIRST TIME RUN TO INITIALIZE THE JUST IN TIME COMPILER
+% Run the ifdiff integration once to compile the code and have a better runtime for a later speed check
+fprintf("Initializing solver %s ... \n", func2str(intIfdiff));
+datahandle = prepareDatahandleForIntegration('pprhs', 'solver', intIfdiff, 'options', intOptions);
+solveODE(datahandle, tspan, x0, p);
+disp("Compilation done...");
 
 %% Solve with IFDIFF
 datahandle = prepareDatahandleForIntegration(rhs, 'solver', intIfdiff, 'options', intOptions);
@@ -53,7 +59,7 @@
 configOld = makeConfig(configNew);
 try
     solIfdiff = solveODE(datahandle, tspan, x0, p);
-catch ME
+catch MEa
     makeConfig(configOld);
     rethrow(ME);
 end
@@ -73,12 +79,46 @@
 
 
 %% Plot solution
+windowTitleX = "3D Solution Trajectories";
 plot_x_ifdiff = deval(solIfdiff, plot_t);
-axSol = plotSol3d([], plot_x_ifdiff, nameIfdiff, lwIfdiff, colorIfdiff, lsIfdiff);
+axSol = plotSol3d([], plot_x_ifdiff, nameIfdiff, lwIfdiff, colorIfdiff, lsIfdiff, windowTitleX);
 
 plot_x_euler = interp1(solEuler.x, solEuler.y', plot_t)';
-plotSol3d(axSol, plot_x_euler, nameEuler, lwEuler, colorEuler, lsEuler);
+plotSol3d(axSol, plot_x_euler, nameEuler, lwEuler, colorEuler, lsEuler, windowTitleX);
 
+%% Plot difference between Euler and IFDIFF solutions
+plot_diff = abs(plot_x_euler - plot_x_ifdiff);
+
+% Plot only every nth entry
+idx = 1:10:length(plot_t);
+
+figure('Name', "Distances of solution states over time");
+for idx_y = 1:3
+    semilogy(plot_t(idx), plot_diff(idx_y, idx), 'LineWidth', 1, 'DisplayName', sprintf('Species %d', idx_y));
+    hold on;
+end
+hold off;
+
+grid on;
+xlabel('Time');
+ylabel('Euler - IFDIFF');
+legend('location', 'northeast');
+title('Difference between Euler and IFDIFF solutions');
+
+%% Plot Euclidean distance between Euler and IFDIFF solutions
+euclidean_diff = vecnorm(plot_x_euler - plot_x_ifdiff, 2, 1);
+
+% Plot only every nth entry
+idx = 1:10:length(plot_t);
+
+figure('Name', "Euclidian distance of solution trajectories over time");
+semilogy(plot_t(idx), euclidean_diff(idx), 'LineWidth', 1, 'DisplayName', 'Euclidean distance');
+
+grid on;
+xlabel('Time');
+ylabel('Euclidean distance');
+legend('location', 'northeast');
+title('Euclidean distance between Euler and IFDIFF solutions');
 
 %% Compute Sensitivity with IFDIFF
 FDstep = generateFDstep(numel(x0), numel(p), 'hy', 1e-6, 'hp', 1e-6, 'ht', 1e-6);
@@ -104,9 +144,25 @@
 
 
 %% Plot Sensitivity
+windowTitle = "Sensitivity and failiure modes of Euler nr. ";
 for idx_y=1:3
-ax = plotSens([], plot_t, plot_sens_ifdiff(idx_y, :), 'sensIFDIFF', lwIfdiff, colorIfdiff, lsIfdiff);
-plotSens(ax, plot_t, plot_sens_euler(idx_y, :), 'sensEuler', lwEuler, colorEuler, lsEuler);
+    windowTitleX = windowTitle + idx_y;
+    % Preferred zoom for each sensitivity
+    switch idx_y
+        case 1
+            zoomX = [79 81];
+            zoomY = [0.32 0.385];
+
+        case 2
+            zoomX = [61.2 61.3];
+            zoomY = [17 30];
+
+        case 3
+            zoomX = [47.5 47.8];
+            zoomY = [9.1 9.7];
+    end
+    ax = plotSens([], plot_t, plot_sens_ifdiff(idx_y, :), 'sensIFDIFF', lwIfdiff, colorIfdiff, lsIfdiff, windowTitleX, zoomX, zoomY);
+    plotSens(ax, plot_t, plot_sens_euler(idx_y, :), 'sensEuler', lwEuler, colorEuler, lsEuler, windowTitleX, zoomX, zoomY);
 end
 
 
@@ -125,14 +181,15 @@
 
 %% Helpers
 function sol_euler = loadEulerOrCompute(fname, vname, rhs, tspan, x0, p, step)
+saveEvery = 100000;
 if isfile(fname)
     tmp = load(fname, vname);
     sol_euler = tmp.(vname);
     fprintf('Loading %s from file %s\n', vname, fname);
     return
 end
-fprintf('Integrating with explicit Euler (might take a while) ...\n')
-sol_euler = explEuler(rhs, tspan, x0, p, step);
+fprintf('Integrating with explicit Euler (might take a while) ...\n');
+sol_euler = explEuler(rhs, tspan, x0, p, step, saveEvery);
 fprintf('Saving result to %s for later reuse.\n', fname);
 save(fname, vname);
 end
@@ -143,7 +200,7 @@ function plotSwitches(ax, switches)
 
 function ax = plotSwitchingFuncOrAlpha(ax, t, x, name, lw, color, yyax, ylab)
 if isempty(ax)
-    f = figure;
+    f = figure('Name', 'Alpha and Switching Function');
     ax = axes(f);
 end
 yyaxis(ax, yyax);
@@ -160,23 +217,41 @@ function plotSwitches(ax, switches)
 ylim(ax, yl);
 end
 
-function ax = plotSens(ax, t, x, name, lw, color, ls)
+function ax = plotSens(ax, t, x, name, lw, color, ls, windowTitleX, zoomX, zoomY)
 if isempty(ax)
-    f = figure;
-    ax = axes(f);
+    f = figure('Name', windowTitleX, 'NumberTitle', 'off');
+    tl = tiledlayout(f, 2, 1);
+    ax = nexttile(tl, 1);
+    axZoom = nexttile(tl, 2);
+    ax.UserData = axZoom; % Store reference to zoom axes
+else
+    axZoom = ax.UserData;
 end
+
+% Full plot
 hold(ax, 'on');
 plot(ax, t, x, 'DisplayName', name, 'LineWidth', lw, 'Color', color, 'LineStyle', ls);
 hold(ax, 'off');
 grid(ax, 'on');
-xlabel(ax, 'Time')
-ylabel(ax, 'Sensitivity')
+xlabel(ax, 'Time');
+ylabel(ax, 'Sensitivity');
 legend(ax, 'location', 'northeast');
+
+% Zoomed plot
+hold(axZoom, 'on');
+plot(axZoom, t, x, 'DisplayName', name, 'LineWidth', lw-1, 'Color', color);
+hold(axZoom, 'off');
+grid(axZoom, 'on');
+xlabel(axZoom, 'Time');
+ylabel(axZoom, 'Sensitivity');
+
+xlim(axZoom, zoomX);
+ylim(axZoom, zoomY);
 end
 
-function ax = plotSol3d(ax, x, name, lw, color, ls)
+function ax = plotSol3d(ax, x, name, lw, color, ls, windowTitleX)
 if isempty(ax)
-    f = figure;
+    f = figure('Name', windowTitleX);
     ax = axes(f);
 end
 hold(ax, 'on');
@@ -192,27 +267,34 @@ function plotSwitches(ax, switches)
 legend(ax, 'location', 'northeast');
 end
 
-function sol = explEuler(rhs, tspan, x0, p, stepsize)
+function sol = explEuler(rhs, tspan, x0, p, stepsize, saveEvery)
 xdim = length(x0);
-stepcount = (tspan(end)-tspan(1)) / stepsize;
-sfac = 0.001; % store factor
-n_out = ceil(stepcount*sfac) + 1;
-
+stepcount = round((tspan(end) - tspan(1)) / stepsize);
+n_out = ceil(stepcount / saveEvery) + 1;
 Xi = reshape(x0, [], 1);
+% Output data structures
 X = zeros(xdim, n_out);
-X(:,1) = Xi;
-
-k = 2; nextout = ceil(1 / sfac);
-for i=2:stepcount
-    Xi = Xi + stepsize * rhs(i*stepsize, Xi, p);
-    if (i == nextout)
-        X(:,k) = Xi; k = k + 1;
-        nextout = nextout + ceil(1 / sfac);
+T = zeros(1, n_out);
+X(:, 1) = Xi;
+T(1) = tspan(1);
+k = 2;
+for i = 1:stepcount
+    t = tspan(1) + (i - 1) * stepsize;
+    Xi = Xi + stepsize * rhs(t, Xi, p);
+    if mod(i, saveEvery) == 0 || i == stepcount
+        X(:, k) = Xi;
+        T(k) = t + stepsize;
+        k = k + 1;
+    end
+    if ~mod(i, max(1, floor(stepcount / 100)))
+        fprintf('.');
     end
-    if ~mod(i, stepcount / 100), fprintf('.'); end
 end
 fprintf('\n');
-T = linspace(tspan(1), tspan(end), n_out);
+% Remove unused preallocated entries
+X = X(:, 1:k-1);
+T = T(1:k-1);
+
 sol.x = T;
 sol.y = X;
 end

toolbox/internal/tools/explEuler.m

@@ -0,0 +1,22 @@
+function sol = explEuler(rhs, tspan, x0, stepsize)
+   xdim = length(x0);  % get dimension
+   stepcount = (tspan(end)-tspan(1))/stepsize;
+   sfac = 0.001; % store factor
+   X = zeros(xdim, ceil(stepcount*sfac)+1);
+   Xi = reshape(x0, [], 1); 
+   X(:,1) = Xi;
+   k = 2; nextout = ceil(1 / sfac);
+   for i=2:stepcount
+      Xi = Xi + stepsize * rhs(i*stepsize, Xi);
+      if (i == nextout)
+         X(:,k) = Xi; k = k + 1; 
+         nextout = nextout + ceil(1 / sfac);
+      end
+      if ~mod(floor(100*i/stepcount), 10), fprintf('.'); end
+   end
+   fprintf('\n')
+   T = linspace(tspan(1), tspan(end), ceil(stepcount*sfac)+1);
+   sol.x = T;
+   sol.y = X;
+end
+