From a11dbc31bf1b9d09a7ed792f1c9968144edddfb3 Mon Sep 17 00:00:00 2001 From: Lev Gromov Date: Mon, 10 Nov 2025 21:26:54 +0100 Subject: [PATCH 01/14] Add equality comparison for branching signatures --- .../functionGeneration/BranchingSignature.m | 30 +++++++++++++++++++ toolbox/internal/solving/compareSignatures.m | 25 ---------------- 2 files changed, 30 insertions(+), 25 deletions(-) delete mode 100644 toolbox/internal/solving/compareSignatures.m diff --git a/toolbox/internal/functionGeneration/BranchingSignature.m b/toolbox/internal/functionGeneration/BranchingSignature.m index 2eb0e2ad..781be0e7 100644 --- a/toolbox/internal/functionGeneration/BranchingSignature.m +++ b/toolbox/internal/functionGeneration/BranchingSignature.m @@ -121,6 +121,7 @@ end end + %% Getters and Setters % Getters retrieve the cached value for str and hash. function str = get.str(this) str = this.strCached; @@ -150,5 +151,34 @@ this.functionIndex = val; this = this.updateCache(); end + + %% Comparison operators + % Equality operator overloading + function tf = eq(obj1, obj2) + tf = false; + % Objects should be of same class, no conversions. + if ~strcmp(class(obj1), class(obj2)) + return + end + + % Check each property + if ~strcmp(obj1.rhsName, obj2.rhsName) + return + end + if ~isequal(obj1.switchCond, obj2.switchCond) + return + end + if ~isequal(obj1.ctrlifIndex, obj2.ctrlifIndex) + return + end + if ~isequal(obj1.functionIndex, obj2.functionIndex) + return + end + tf = true; + end + + function tf = ne(obj1, obj2) + tf = ~eq(obj1, obj2); + end end end diff --git a/toolbox/internal/solving/compareSignatures.m b/toolbox/internal/solving/compareSignatures.m deleted file mode 100644 index 338b8205..00000000 --- a/toolbox/internal/solving/compareSignatures.m +++ /dev/null @@ -1,25 +0,0 @@ -function flag = compareSignatures(ctrlif_index1, ctrlif_index2, function_index1, function_index2, switch_cond1, switch_cond2) -% flag = compareSignatures(ctrlif_index1, ctrlif_index2, function_index1, function_index2, switch_cond1, switch_cond2) -% -% Compares two signatures. Returns true if they are equal, false if not. - -flag = true; - -if ~all(switch_cond1 == switch_cond2) || ~all(ctrlif_index1 == ctrlif_index2) - flag = false; - return -end - -if length(function_index1) ~= length(function_index2) - flag = false; - return -else - for i=1:length(function_index1) - if ~all(function_index1{i}==function_index2{i}) - flag = false; - break - end - end -end - -end \ No newline at end of file From c9caaf42ffc9a49de3424902d6fb3efa492b4b16 Mon Sep 17 00:00:00 2001 From: Lev Gromov Date: Mon, 10 Nov 2025 21:43:56 +0100 Subject: [PATCH 02/14] Store signatures as BranchingSignature objects in datahandle Also store Filippov signatures as a cell consisting of the two signatures belonging to the overlapping models. --- .../internal/sensitivities/getModelFunction.m | 11 +- .../solving/chatteringGetSignatures.m | 88 +++---- toolbox/internal/solving/ctrlif.m | 242 +++++++++--------- .../initDatahandleFields_SWP_detection.m | 5 +- .../solving/slidingFilippovRHS_oneSwitch.m | 14 +- .../solveODE_cutSteps_solution_until_t2.m | 24 +- toolbox/internal/solving/solveODE_firstTime.m | 52 ++-- .../solving/solveODE_prepareNextStage.m | 5 +- .../solving/solveODE_setFilippovRHS.m | 22 +- 9 files changed, 223 insertions(+), 240 deletions(-) diff --git a/toolbox/internal/sensitivities/getModelFunction.m b/toolbox/internal/sensitivities/getModelFunction.m index 8233ccfd..44d803b5 100644 --- a/toolbox/internal/sensitivities/getModelFunction.m +++ b/toolbox/internal/sensitivities/getModelFunction.m @@ -17,16 +17,7 @@ data = datahandle.getData(); factory = data.codeGen.modelFactory; -rhsName = data.mtreeplus{2, 1}; -dataSignature = data.SWP_detection.signature; -switchCond = dataSignature.switch_cond{modelStage}; -ctrlifIndex = dataSignature.ctrlif_index{modelStage}; -functionIndex = dataSignature.function_index{modelStage}; -signature = BranchingSignature( ... - rhsName, ... - switchCond, ... - ctrlifIndex, ... - functionIndex); +signature = data.SWP_detection.signature{modelStage}; [functionHandle, collisionIndex] = factory.findExisting(signature); if isempty(functionHandle) diff --git a/toolbox/internal/solving/chatteringGetSignatures.m b/toolbox/internal/solving/chatteringGetSignatures.m index be78a17c..709d74a1 100644 --- a/toolbox/internal/solving/chatteringGetSignatures.m +++ b/toolbox/internal/solving/chatteringGetSignatures.m @@ -1,22 +1,22 @@ function [signature1, signature2] = chatteringGetSignatures(datahandle, t) % [signature1, signature2] = chatteringGetSignatures(datahandle, t) -% -% For a chattering solution, determines the signatures that appeared +% +% For a chattering solution, determines the signatures that appeared % since timepoint t. If there are more than two, throws an error since the % program is not capable of computing a Filippov solution in that case. % % % INPUT: -% datahandle: datahandle containing switching point information +% datahandle: datahandle containing switching point information % handle -% +% % t: timepoint, ignore signatures before this point % scalar % % OUTPUT: % signature1: Signature that has been encountered. % struct -% +% % signature2: Signature that has been encountered. % struct @@ -27,53 +27,55 @@ % Assume switching points are sorted. If not, something went very wrong... idxChatterStart = find(switchingpoints >= t, 1); -signatures = data.SWP_detection.signature; -ctrlif_indices = signatures.ctrlif_index(idxChatterStart:end); -function_indices = signatures.function_index(idxChatterStart:end); -switch_conds = signatures.switch_cond(idxChatterStart:end); - -k = length(ctrlif_indices); - -uniqueChatteringSignatures.ctrlif_index = ctrlif_indices(1); -uniqueChatteringSignatures.function_index = function_indices(1); -uniqueChatteringSignatures.switch_cond = switch_conds(1); - -for i=2:k - exists = false; - ctrlif_index = ctrlif_indices{i}; - function_index = function_indices{i}; - switch_cond = switch_conds{i}; +signatures = data.SWP_detection.signature(idxChatterStart:end); +% Make sure we don't consider any signatures from previous Filippov modes. +% Filippov signatures are stored as cells with two signatures. +isFilippov = cellfun(@iscell, signatures); +if any(isFilippov) + % Unsure if this case can be handled gracefully, for now abort. + throw(filippovSignatureChatteringException); +end - for j=1:length(uniqueChatteringSignatures.ctrlif_index) - ctrlif_index_ref = uniqueChatteringSignatures.ctrlif_index{j}; - function_index_ref = uniqueChatteringSignatures.function_index{j}; - switch_cond_ref = uniqueChatteringSignatures.switch_cond{j}; - if compareSignatures( ... - ctrlif_index, ctrlif_index_ref, ... - function_index, function_index_ref, ... - switch_cond, switch_cond_ref) - exists = true; +% Find unique signatures among chattering. Stop and throw if we find more than two. +uniqueChatteringSignatures = signatures(1); +for idxNew=2:length(signatures) + isNewSignature = true; + for idxOld=1:length(uniqueChatteringSignatures) + if signatures{idxNew} == signatures{idxOld} + % Signature already included. + isNewSignature = false; + break end end - % signature not seen yet - if ~exists - uniqueChatteringSignatures.ctrlif_index{end+1} = ctrlif_index; - uniqueChatteringSignatures.function_index{end+1} = function_index; - uniqueChatteringSignatures.switch_cond{end+1} = switch_cond; + + if isNewSignature + % New signature, add to list. + uniqueChatteringSignatures{end + 1} = signatures{idxNew}; %#ok + end + + if length(uniqueChatteringSignatures) > 2 + break end end -if length(uniqueChatteringSignatures.ctrlif_index) ~= 2 - error('IFDIFF:ChatteringInvolvesSeveralSwitches', 'Chattering does not involve exactly two signatures.\n'); + +if length(uniqueChatteringSignatures) ~= 2 + throw(invalidNumberOfChatteringSwitches); end % output -signature1.ctrlif_index = uniqueChatteringSignatures.ctrlif_index{1}; -signature1.function_index = uniqueChatteringSignatures.function_index{1}; -signature1.switch_cond = uniqueChatteringSignatures.switch_cond{1}; +[signature1, signature2] = uniqueChatteringSignatures{:}; +end -signature2.ctrlif_index = uniqueChatteringSignatures.ctrlif_index{2}; -signature2.function_index = uniqueChatteringSignatures.function_index{2}; -signature2.switch_cond = uniqueChatteringSignatures.switch_cond{2}; +%% Exceptions +function e = filippovSignatureChatteringException +msg = [ ... + 'Unable to determine chattering signature:\n' ... + 'One of the chattering signature candidates is a Filippov signature.\n']; +e = MException('IFDIFF:ChatteringInvolvesFilippov', msg); +end +function e = invalidNumberOfChatteringSwitches +msg = 'Chattering does not involve exactly two signatures.\n'; +e = MException('IFDIFF:ChatteringInvolvesSeveralSwitches', msg); end diff --git a/toolbox/internal/solving/ctrlif.m b/toolbox/internal/solving/ctrlif.m index 98e3f8d9..7d36a81c 100644 --- a/toolbox/internal/solving/ctrlif.m +++ b/toolbox/internal/solving/ctrlif.m @@ -1,116 +1,128 @@ function y = ctrlif(num, truepart, falsepart, ctrlif_index, function_index, datahandle) - % function for: - % - % - supervising signature during integration (.active_SWP_detection = 1) - % - forced evaluation condition during extending of the ODE solution object until switch - % - receiving the siganture at any timepoint t (.getSignature = 1) - % - forced evaluation with a given sigature for sensitivity calculation - % - % by the ctrilf one gets the condition, ctrlif_index, function_index - % those vectors are preallocated using .ctrlifCounter - - config = makeConfig(); - data = datahandle.getData(); - - condition = (num >= 0); - - switch data.caseCtrlif - case config.caseCtrlif.default - if condition - y = truepart; - else - y = falsepart; - end - - case config.caseCtrlif.forcedBranching - % active_SWP_detection: during integration used in extendODE, solveODE - % forced branching method: always return the initial true/false value - - data.forcedBranching.ctrlifCounter = data.forcedBranching.ctrlifCounter + 1; - - data.forcedBranching.switch_cond(data.forcedBranching.ctrlifCounter) = condition; - data.forcedBranching.ctrlif_index(data.forcedBranching.ctrlifCounter) = ctrlif_index; - data.forcedBranching.function_index{data.forcedBranching.ctrlifCounter} = function_index; - - - forced_evaluation_cond = data.forcedBranching.switch_cond_forcedBranching(data.forcedBranching.ctrlifCounter); - - if forced_evaluation_cond - y = truepart; - else - y = falsepart; - end - - % ctrlifCounter for preallocation of signature - - if length(data.forcedBranching.switch_cond_forcedBranching) == data.forcedBranching.ctrlifCounter - data.forcedBranching.ctrlifCounter = 0; - end - - case config.caseCtrlif.extendODEuntilSwitch - - data.forcedBranching.ctrlifCounter = data.forcedBranching.ctrlifCounter + 1; - - forced_evaluation_cond = data.forcedBranching.switch_cond_forcedBranching(data.forcedBranching.ctrlifCounter); - - if forced_evaluation_cond - y = truepart; - else - y = falsepart; - end - - % ctrlifCounter for preallocation of signature - - if length(data.forcedBranching.switch_cond_forcedBranching) == data.forcedBranching.ctrlifCounter - data.forcedBranching.ctrlifCounter = 0; - end - - case config.caseCtrlif.getSignature - data.getSignature.ctrlifCounter = data.getSignature.ctrlifCounter + 1; - - data.getSignature.switch_cond = ctrlif_getSignaturePreallocation(condition, data.getSignature.switch_cond, data.getSignature.ctrlifCounter, 1); - data.getSignature.ctrlif_index = ctrlif_getSignaturePreallocation(ctrlif_index, data.getSignature.ctrlif_index, data.getSignature.ctrlifCounter, 1); - data.getSignature.function_index = ctrlif_getSignaturePreallocation(function_index, data.getSignature.function_index, data.getSignature.ctrlifCounter, 2); - - if condition - y = truepart; - else - y = falsepart; - end - case config.caseCtrlif.getSignatureChange - data.getSignature.ctrlifCounter = data.getSignature.ctrlifCounter + 1; - - data.getSignature.switch_cond = ctrlif_getSignaturePreallocation(condition, data.getSignature.switch_cond, data.getSignature.ctrlifCounter, 1); - data.getSignature.ctrlif_index = ctrlif_getSignaturePreallocation(ctrlif_index, data.getSignature.ctrlif_index, data.getSignature.ctrlifCounter, 1); - data.getSignature.function_index = ctrlif_getSignaturePreallocation(function_index, data.getSignature.function_index, data.getSignature.ctrlifCounter, 2); - - forced_evaluation_cond = data.getSignature.switch_cond_forcedBranching(data.getSignature.ctrlifCounter); - if forced_evaluation_cond - y = truepart; - else - y = falsepart; - end - case config.caseCtrlif.computeSensitivities - - % evaluate RHS with a given signature for calculating the sensitivities - - modelStage = data.computeSensitivity.modelStage; - - data.computeSensitivity.ctrlifCounter = data.computeSensitivity.ctrlifCounter + 1; - signature = data.SWP_detection.signature.switch_cond{modelStage}; - forced_evaluation_cond = signature(data.computeSensitivity.ctrlifCounter); - - if forced_evaluation_cond - y = truepart; - else - y = falsepart; - end - - % ctrlifCounter for preallocation of signature - if length(signature) == data.computeSensitivity.ctrlifCounter - data.computeSensitivity.ctrlifCounter = 0; - end - end - - datahandle.setData(data); -end \ No newline at end of file +% function for: +% +% - supervising signature during integration (.active_SWP_detection = 1) +% - forced evaluation condition during extending of the ODE solution object until switch +% - receiving the siganture at any timepoint t (.getSignature = 1) +% - forced evaluation with a given sigature for sensitivity calculation +% +% by the ctrilf one gets the condition, ctrlif_index, function_index +% those vectors are preallocated using .ctrlifCounter + +config = makeConfig(); +data = datahandle.getData(); + +condition = (num >= 0); + +switch data.caseCtrlif + case config.caseCtrlif.default + if condition + y = truepart; + else + y = falsepart; + end + + case config.caseCtrlif.forcedBranching + % active_SWP_detection: during integration used in extendODE, solveODE + % forced branching method: always return the initial true/false value + + data.forcedBranching.ctrlifCounter = data.forcedBranching.ctrlifCounter + 1; + + data.forcedBranching.switch_cond(data.forcedBranching.ctrlifCounter) = condition; + data.forcedBranching.ctrlif_index(data.forcedBranching.ctrlifCounter) = ctrlif_index; + data.forcedBranching.function_index{data.forcedBranching.ctrlifCounter} = function_index; + + + forced_evaluation_cond = data.forcedBranching.switch_cond_forcedBranching(data.forcedBranching.ctrlifCounter); + + if forced_evaluation_cond + y = truepart; + else + y = falsepart; + end + + % ctrlifCounter for preallocation of signature + + if length(data.forcedBranching.switch_cond_forcedBranching) == data.forcedBranching.ctrlifCounter + data.forcedBranching.ctrlifCounter = 0; + end + + case config.caseCtrlif.extendODEuntilSwitch + + data.forcedBranching.ctrlifCounter = data.forcedBranching.ctrlifCounter + 1; + + forced_evaluation_cond = data.forcedBranching.switch_cond_forcedBranching(data.forcedBranching.ctrlifCounter); + + if forced_evaluation_cond + y = truepart; + else + y = falsepart; + end + + % ctrlifCounter for preallocation of signature + + if length(data.forcedBranching.switch_cond_forcedBranching) == data.forcedBranching.ctrlifCounter + data.forcedBranching.ctrlifCounter = 0; + end + + case config.caseCtrlif.getSignature + data.getSignature.ctrlifCounter = data.getSignature.ctrlifCounter + 1; + + data.getSignature.switch_cond = ctrlif_getSignaturePreallocation(condition, data.getSignature.switch_cond, data.getSignature.ctrlifCounter, 1); + data.getSignature.ctrlif_index = ctrlif_getSignaturePreallocation(ctrlif_index, data.getSignature.ctrlif_index, data.getSignature.ctrlifCounter, 1); + data.getSignature.function_index = ctrlif_getSignaturePreallocation(function_index, data.getSignature.function_index, data.getSignature.ctrlifCounter, 2); + + if condition + y = truepart; + else + y = falsepart; + end + case config.caseCtrlif.getSignatureChange + data.getSignature.ctrlifCounter = data.getSignature.ctrlifCounter + 1; + + data.getSignature.switch_cond = ctrlif_getSignaturePreallocation(condition, data.getSignature.switch_cond, data.getSignature.ctrlifCounter, 1); + data.getSignature.ctrlif_index = ctrlif_getSignaturePreallocation(ctrlif_index, data.getSignature.ctrlif_index, data.getSignature.ctrlifCounter, 1); + data.getSignature.function_index = ctrlif_getSignaturePreallocation(function_index, data.getSignature.function_index, data.getSignature.ctrlifCounter, 2); + + forced_evaluation_cond = data.getSignature.switch_cond_forcedBranching(data.getSignature.ctrlifCounter); + if forced_evaluation_cond + y = truepart; + else + y = falsepart; + end + case config.caseCtrlif.computeSensitivities + + % evaluate RHS with a given signature for calculating the sensitivities + + modelStage = data.computeSensitivity.modelStage; + + data.computeSensitivity.ctrlifCounter = data.computeSensitivity.ctrlifCounter + 1; + signature = data.SWP_detection.signature{modelStage}; + + % Cannot evaluate model with Filippov signature (i.e. two signatures). + if iscell(signature) + throw(evaluationWithFilippovSignatureException) + end + + forced_evaluation_cond = signature.switchCond(data.computeSensitivity.ctrlifCounter); + + if forced_evaluation_cond + y = truepart; + else + y = falsepart; + end + + % ctrlifCounter for preallocation of signature + if length(signature.switchCond) == data.computeSensitivity.ctrlifCounter + data.computeSensitivity.ctrlifCounter = 0; + end +end + +datahandle.setData(data); +end + +%% Exceptions +function e = evaluationWithFilippovSignatureException +msg = 'Cannot evaluate model in computeSensitivity mode for Filippov signature.\n'; +e = MException('IFDIFF:Ctrlif:ModelEvaluationForFilippovSignature', msg); +end diff --git a/toolbox/internal/solving/initDatahandleFields_SWP_detection.m b/toolbox/internal/solving/initDatahandleFields_SWP_detection.m index 51ed2463..da529cc9 100644 --- a/toolbox/internal/solving/initDatahandleFields_SWP_detection.m +++ b/toolbox/internal/solving/initDatahandleFields_SWP_detection.m @@ -26,13 +26,10 @@ SWP_detection.function_index_t3 = []; SWP_detection.signature = {}; - SWP_detection.signature.function_index = {}; - SWP_detection.signature.ctrlif_index = {}; - SWP_detection.signature.switch_cond = {}; SWP_detection.solution_until_t1 = {}; SWP_detection.solution_until_t2 = {}; SWP_detection.solution_until_t3 = {}; SWP_detection.uniqueSwEnumeration = 0; -end \ No newline at end of file +end diff --git a/toolbox/internal/solving/slidingFilippovRHS_oneSwitch.m b/toolbox/internal/solving/slidingFilippovRHS_oneSwitch.m index 4ce0007b..7b95a7fc 100644 --- a/toolbox/internal/solving/slidingFilippovRHS_oneSwitch.m +++ b/toolbox/internal/solving/slidingFilippovRHS_oneSwitch.m @@ -38,15 +38,15 @@ % Heidelberg, 2020. % evaluate RHS for signature where chattering switchingFunction is positive -ctrlif_index_fminus = signature_fminus.ctrlif_index; -function_index_fminus = signature_fminus.function_index; -switch_cond_fminus = signature_fminus.switch_cond; +ctrlif_index_fminus = signature_fminus.ctrlifIndex; +function_index_fminus = signature_fminus.functionIndex; +switch_cond_fminus = signature_fminus.switchCond; f_minus = evaluateRHSForSignature(datahandle, t, y, p, ctrlif_index_fminus, function_index_fminus, switch_cond_fminus); % evaluate RHS for signature where chattering switchingFunction is negative -ctrlif_index_fplus = signature_fplus.ctrlif_index; -function_index_fplus = signature_fplus.function_index; -switch_cond_fplus = signature_fplus.switch_cond; +ctrlif_index_fplus = signature_fplus.ctrlifIndex; +function_index_fplus = signature_fplus.functionIndex; +switch_cond_fplus = signature_fplus.switchCond; f_plus = evaluateRHSForSignature(datahandle, t, y, p, ctrlif_index_fplus, function_index_fplus, switch_cond_fplus); @@ -80,4 +80,4 @@ dy = alpha*f_minus + (1-alpha)*f_plus; -end \ No newline at end of file +end diff --git a/toolbox/internal/solving/solveODE_cutSteps_solution_until_t2.m b/toolbox/internal/solving/solveODE_cutSteps_solution_until_t2.m index ed39fcf2..af05b5df 100644 --- a/toolbox/internal/solving/solveODE_cutSteps_solution_until_t2.m +++ b/toolbox/internal/solving/solveODE_cutSteps_solution_until_t2.m @@ -1,7 +1,7 @@ function solveODE_cutSteps_solution_until_t2(datahandle, t_cut) -% Removes integration step that exceed time 't_cut' from ode sol object in +% Removes integration step that exceed time 't_cut' from ode sol object in % datahandle. -% +% % INPUT: % 'datahandle': datahandle containing the integration and switching data. % handle @@ -19,7 +19,7 @@ function solveODE_cutSteps_solution_until_t2(datahandle, t_cut) % % Code adapted from solveODE_solution_until_t1.m -data = datahandle.getData(); +data = datahandle.getData(); % Step 1: In solution object, cut steps past time t_cut cutSteps_solution(datahandle, 'solution_until_t2', t_cut); @@ -28,16 +28,18 @@ function solveODE_cutSteps_solution_until_t2(datahandle, t_cut) SWP_detection = data.SWP_detection; switchingpoints = cell2mat(SWP_detection.switchingpoints); -SWP_detection.switchingpoints = SWP_detection.switchingpoints(switchingpoints<=t_cut); -SWP_detection.switchingFunction = SWP_detection.switchingFunction(switchingpoints<=t_cut); -SWP_detection.jumpFunction = SWP_detection.jumpFunction(switchingpoints<=t_cut); +indicesCut = switchingpoints <= t_cut; + +SWP_detection.switchingpoints = SWP_detection.switchingpoints(indicesCut); +SWP_detection.switchingFunction = SWP_detection.switchingFunction(indicesCut); +SWP_detection.jumpFunction = SWP_detection.jumpFunction(indicesCut); t2 = data.SWP_detection.solution_until_t2.x(end); x2 = data.SWP_detection.solution_until_t2.y(:,end); SWP_detection.t2 = t2; SWP_detection.x2 = {x2, x2}; data.SWP_detection = SWP_detection; -datahandle.setData(data); +datahandle.setData(data); % adjust also signature & co. at t2 [switch_cond_t2, ctrlif_index_t2, function_index_all_t2] = ctrlif_getSignature(... datahandle, ... @@ -48,14 +50,10 @@ function solveODE_cutSteps_solution_until_t2(datahandle, t_cut) data.SWP_detection.function_index_t2 = function_index_all_t2; % note: signature keeps length(switchingpoints)+1 entries, one entry for % the initial signature -data.SWP_detection.signature.ctrlif_index = data.SWP_detection.signature.ctrlif_index([true, switchingpoints<=t_cut]); -data.SWP_detection.signature.function_index = data.SWP_detection.signature.function_index([true, switchingpoints<=t_cut]); -data.SWP_detection.signature.switch_cond = data.SWP_detection.signature.switch_cond([true, switchingpoints<=t_cut]); +data.SWP_detection.signature = data.SWP_detection.signature([true, indicesCut]); datahandle.setData(data); % adjust convexification data cutSteps_convexification(datahandle, t_cut) - - -end \ No newline at end of file +end diff --git a/toolbox/internal/solving/solveODE_firstTime.m b/toolbox/internal/solving/solveODE_firstTime.m index d36b6320..60dbe430 100644 --- a/toolbox/internal/solving/solveODE_firstTime.m +++ b/toolbox/internal/solving/solveODE_firstTime.m @@ -27,7 +27,7 @@ function solveODE_firstTime(datahandle) % get maxStep size by default defaultMaxStep = 0.1*abs(data.SWP_detection.tspan(1)-data.SWP_detection.tspan(2)); -data.integratorSettings.optionsForcedBranching.InitialStep = defaultMaxStep; +data.integratorSettings.optionsForcedBranching.InitialStep = defaultMaxStep; data.caseCtrlif = config.caseCtrlif.forcedBranching; datahandle.setData(data) % initialisation of nIntegrationSteps s.th. the first iteration of while is executed @@ -38,11 +38,11 @@ function solveODE_firstTime(datahandle) % while loop as long as there are only 2 timpoints or less during the integration % algorithm while nIntegrationSteps <= 2 - + % getData since the ode options may have been changed. data = datahandle.getData(); - - + + % caution: .optionsForcedBranching contains analyseSignature(...,datahandle). % switching point detection magic is done in analyseSignature(...,datahandle). z = data.integratorSettings.numericIntegrator(... @@ -50,51 +50,33 @@ function solveODE_firstTime(datahandle) data.SWP_detection.tspan,... data.SWP_detection.initialvalues, ... data.integratorSettings.optionsForcedBranching); - + data = datahandle.getData(); data.SWP_detection.solution_until_t3 = z; - + nIntegrationSteps = length(data.SWP_detection.solution_until_t3.x); - + % reduce max step size if number of timepoints not high enough if nIntegrationSteps <= 2 - + data.integratorSettings.optionsForcedBranching.InitialStep = data.integratorSettings.optionsForcedBranching.InitialStep/2; %warning('First integration step contains switch. odeextend would fail, reduce step size only for first integration step') end - + % save new MaxStep size to datahandle. - datahandle.setData(data); + datahandle.setData(data); end % get signature data.SWP_detection.t1 = data.SWP_detection.solution_until_t3.x(end-1); data.SWP_detection.t3 = data.SWP_detection.solution_until_t3.x(end); -% signature matrix -data.SWP_detection.signature.function_index{1} = data.forcedBranching.function_index_forcedBranching; -data.SWP_detection.signature.ctrlif_index{1} = data.forcedBranching.ctrlif_index_forcedBranching; -data.SWP_detection.signature.switch_cond{1} = data.forcedBranching.switch_cond_forcedBranching; -datahandle.setData(data); - +firstModelSignature = BranchingSignature( ... + data.mtreeplus{2, 1}, ... + data.forcedBranching.switch_cond_forcedBranching, ... + data.forcedBranching.ctrlif_index_forcedBranching, ... + data.forcedBranching.function_index_forcedBranching); +data.SWP_detection.signature{1} = firstModelSignature; +datahandle.setData(data); end % first integraion round done, if no switch occured, .t3 = tspan(2) - - - - - - - - - - - - - - - - - - - diff --git a/toolbox/internal/solving/solveODE_prepareNextStage.m b/toolbox/internal/solving/solveODE_prepareNextStage.m index 15416dbf..e40ded77 100644 --- a/toolbox/internal/solving/solveODE_prepareNextStage.m +++ b/toolbox/internal/solving/solveODE_prepareNextStage.m @@ -16,9 +16,8 @@ function solveODE_prepareNextStage(datahandle) % Get the new signature at t2, without forced branching [switch_cond, ctrlif_index, function_index] = ctrlif_getSignature(datahandle, t, xPlus); -data.SWP_detection.signature.switch_cond{end+1} = switch_cond; -data.SWP_detection.signature.ctrlif_index{end+1} = ctrlif_index; -data.SWP_detection.signature.function_index{end+1} = function_index; +signatureNew = BranchingSignature(data.mtreeplus{2, 1}, switch_cond, ctrlif_index, function_index); +data.SWP_detection.signature{end + 1} = signatureNew; datahandle.setData(data); end diff --git a/toolbox/internal/solving/solveODE_setFilippovRHS.m b/toolbox/internal/solving/solveODE_setFilippovRHS.m index 0ee7526c..1306e33e 100644 --- a/toolbox/internal/solving/solveODE_setFilippovRHS.m +++ b/toolbox/internal/solving/solveODE_setFilippovRHS.m @@ -1,9 +1,9 @@ function solveODE_setFilippovRHS(datahandle) -% Sets datahandle.sliding.filippovRHS to a RHS that allows to slide/run +% Sets datahandle.sliding.filippovRHS to a RHS that allows to slide/run % on the zero-manifold associated to a function that is inconsistently switching. % % INPUT: -% 'datahandle': datahandle containing the integration and +% 'datahandle': datahandle containing the integration and % switching data. % % OUTPUT: @@ -25,7 +25,7 @@ function solveODE_setFilippovRHS(datahandle) % only two signatures involved, so we can assume the last switching ctrlif % to be the chattering one sliding_index = datahandle.getData().SWP_detection.switchingIndices(end); -if signature1.switch_cond(sliding_index) == 1 +if signature1.switchCond(sliding_index) == 1 signature_fplus = signature1; signature_fminus = signature2; else @@ -33,28 +33,30 @@ function solveODE_setFilippovRHS(datahandle) signature_fminus = signature1; end - % cut steps solveODE_cutSteps_solution_until_t2(datahandle, t) +data = datahandle.getData(); + +% Store signature of Filippov model as combined signatures for sensitivities. +data.SWP_detection.signature{end} = {signature_fplus, signature_fminus}; + % get the chattering switch's switching function: -% We ensured there was only a single switch active during chattering, +% We ensured there was only a single switch active during chattering, % so we just pick the switching function of the last switching event. -switchingFunction = datahandle.getData().SWP_detection.switchingfunctionhandles{end}; +switchingFunction = data.SWP_detection.switchingfunctionhandles{end}; filippov_rhs = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch(datahandle, signature_fplus, signature_fminus, switchingFunction,t,y,p); % set filippov rhs and store some other info -data = datahandle.getData(); data.sliding.filippov_rhs = filippov_rhs; data.sliding.index = sliding_index; data.sliding.signature_fplus = signature_fplus; data.sliding.signature_fminus = signature_fminus; + datahandle.setData(data); -% message if config.debugMode fprintf("Entered Filippov regime.\n"); end - -end \ No newline at end of file +end From a4aa57479d6f14dc4db6a1426e106f1012304595 Mon Sep 17 00:00:00 2001 From: Lev Gromov Date: Mon, 10 Nov 2025 21:45:09 +0100 Subject: [PATCH 03/14] Treat exit from Filippov mode as proper switch in datahandle --- toolbox/solveODE.m | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/toolbox/solveODE.m b/toolbox/solveODE.m index bfa952cc..7bba4a6b 100644 --- a/toolbox/solveODE.m +++ b/toolbox/solveODE.m @@ -60,12 +60,24 @@ if checkForSwitchingIndices(datahandle) error('IFDIFF:switchInFilippov', 'Switching during or right after Filippov mode.'); end - if datahandle.getData().SWP_detection.t2 == tspan(2) + data = datahandle.getData(); + if data.SWP_detection.t2 == tspan(2) break; % reached end of time span else if config.debugMode fprintf('Left Filippov regime.\n'); end + % Left the Filippov model and continue integration with new model: + % Need to set switching point, switching function and new signature. + % TODO: Do we need an accurate exit time of the Filippov model, e.g. for sensitivities? + % TODO: Same goes for the switching function... + % For now, just take the time where we ended up leaving and same switching function as in entry into Filippov. + data = datahandle.getData(); + data.SWP_detection.switchingpoints{end + 1} = data.SWP_detection.t2; + % TODO: This in particular seems very wrong! + data.SWP_detection.switchingFunction{end + 1} = data.SWP_detection.switchingFunction{end}; + data.SWP_detection.jumpFunction{end + 1} = []; + datahandle.setData(data); solveODE_prepareNextStage(datahandle); end end From f5f70576af791f1f31c5a5244b7f7614506bcfd2 Mon Sep 17 00:00:00 2001 From: Lev Gromov Date: Mon, 10 Nov 2025 23:19:34 +0100 Subject: [PATCH 04/14] Fix RHS eval for signature not resetting the ctrlif mode --- toolbox/internal/solving/evaluateRHSForSignature.m | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/toolbox/internal/solving/evaluateRHSForSignature.m b/toolbox/internal/solving/evaluateRHSForSignature.m index fa3dc9f0..d96ec6c2 100644 --- a/toolbox/internal/solving/evaluateRHSForSignature.m +++ b/toolbox/internal/solving/evaluateRHSForSignature.m @@ -39,6 +39,7 @@ % counter from the ctrlif data.forcedBranching.ctrlifCounter = 0; +oldCaseCtrlif = data.caseCtrlif; data.caseCtrlif = config.caseCtrlif.forcedBranching; % case forced branching datahandle.setData(data); @@ -47,5 +48,6 @@ % cleanup data = datahandle.getData(); -data.caseCtrlif = config.caseCtrlif.default; % default case -end \ No newline at end of file +data.caseCtrlif = oldCaseCtrlif; +datahandle.setData(data); +end From 00158a57a929b758b22dc64caa9a378dc7f3fd96 Mon Sep 17 00:00:00 2001 From: Lev Gromov Date: Mon, 10 Nov 2025 23:26:16 +0100 Subject: [PATCH 05/14] WIP: Add support for Filippov sensitivities Currently only works for systems where the Filippov sliding mode is entered, but never exited. --- .../internal/sensitivities/getGy_Gp_update.m | 173 ++++++++++-------- .../sensitivities/solveDisturbed_Gp.m | 49 +++-- .../sensitivities/solveDisturbed_Gy.m | 43 +++-- toolbox/internal/sensitivities/solveVDE_Gp.m | 39 ++-- toolbox/internal/sensitivities/solveVDE_Gy.m | 36 ++-- 5 files changed, 203 insertions(+), 137 deletions(-) diff --git a/toolbox/internal/sensitivities/getGy_Gp_update.m b/toolbox/internal/sensitivities/getGy_Gp_update.m index 838ed61b..4005f5dc 100644 --- a/toolbox/internal/sensitivities/getGy_Gp_update.m +++ b/toolbox/internal/sensitivities/getGy_Gp_update.m @@ -1,91 +1,106 @@ function Updates = getGy_Gp_update(datahandle, startModel, endModel, Gp_flag, options) - % Updates = getGy_Gp_update(datahandle, startModel, endModel, Gp_flag, options) - % - % Calculates the update-matrices for the sensitivitiy calculation with END_piecewise or VDE. - % - % INPUT: datahandle - datahandle you get from the integration process with solveODE - % amountG - amount of intermediate G-matrices that already have been calculated and saved from previous function calls. - % (If amountG = 1, it means that the intermediates were only initialised and all the matrices have still to be calculated.) - % modelNum - number of the model until which the update-matrices need to be calculated - % Gp_flag - flag that is true if the sensitivities with respect to the parameters should be calculated - % options - struct that has the informations for the stepsizes for END, the integrator and the method - % - % OUTPUT: Updates - struct that has he update-matrices needed for the sensitivities w.r.t. y0 and optionally w.r.t. p - - % Initialize - data = datahandle.getData(); - parameters = data.SWP_detection.parameters; - switches = data.computeSensitivity.switches_extended; - switches_left = data.computeSensitivity.switches_extended_left; - y_to_switches = data.computeSensitivity.y_to_switches; - y_to_switches_left = data.computeSensitivity.y_to_switches_left; - switching_functions = data.SWP_detection.switchingFunction; - jump_functions = data.SWP_detection.jumpFunction; - functionRHS = data.integratorSettings.preprocessed_rhs; - - FDstep = options.FDstep; - - dim_y = data.computeSensitivity.dim_y; - dim_p = data.computeSensitivity.dim_p; - unit = eye(dim_y); - - numModels = endModel - startModel; - Updates.Uy_new = cell(1, numModels); +% Updates = getGy_Gp_update(datahandle, startModel, endModel, Gp_flag, options) +% +% Calculates the update-matrices for the sensitivitiy calculation with END_piecewise or VDE. +% +% INPUT: datahandle - datahandle you get from the integration process with solveODE +% amountG - amount of intermediate G-matrices that already have been calculated and saved from previous function calls. +% (If amountG = 1, it means that the intermediates were only initialised and all the matrices have still to be calculated.) +% modelNum - number of the model until which the update-matrices need to be calculated +% Gp_flag - flag that is true if the sensitivities with respect to the parameters should be calculated +% options - struct that has the informations for the stepsizes for END, the integrator and the method +% +% OUTPUT: Updates - struct that has he update-matrices needed for the sensitivities w.r.t. y0 and optionally w.r.t. p - if Gp_flag - Updates.Up_new = cell(1, numModels); - h_p = fdStep_getH_p(FDstep, parameters); +% Initialize +data = datahandle.getData(); +parameters = data.SWP_detection.parameters; +switches = data.computeSensitivity.switches_extended; +switches_left = data.computeSensitivity.switches_extended_left; +y_to_switches = data.computeSensitivity.y_to_switches; +y_to_switches_left = data.computeSensitivity.y_to_switches_left; +switching_functions = data.SWP_detection.switchingFunction; +jump_functions = data.SWP_detection.jumpFunction; +functionRHS = data.integratorSettings.preprocessed_rhs; + +FDstep = options.FDstep; + +dim_y = data.computeSensitivity.dim_y; +dim_p = data.computeSensitivity.dim_p; +unit = eye(dim_y); + +numModels = endModel - startModel; +Updates.Uy_new = cell(1, numModels); + +if Gp_flag + Updates.Up_new = cell(1, numModels); + h_p = fdStep_getH_p(FDstep, parameters); +end + +% Fix the model and set the model number for which model you want to evaluate the RHS +data.computeSensitivity.modelStage = startModel; +datahandle.setData(data); + +% Calculate the remaining update matrices for the update formula until modelNum. +for i = startModel : (endModel-1) + switchingFunction = switching_functions{i}; + jumpFunction = jump_functions{i}; + ts = switches(i+1); + tsMinus = switches_left(i+1); + yMinus = y_to_switches_left(:, i+1); + yPlus = y_to_switches(:, i+1); + + h_y = fdStep_getH_y(FDstep, yMinus); + h_t = fdStep_getH_t(FDstep, tsMinus); + + % Calculate the derivatives of the switching functions w.r.t. y, t (and p if necessary) + del_sigmay = del_f_del_y(datahandle, switchingFunction, tsMinus, yMinus, parameters, h_y); + del_sigmat = del_f_del_t(datahandle, switchingFunction, tsMinus, yMinus, parameters, h_t); + diff_sigmat = del_sigmat + del_sigmay * functionRHS(datahandle, tsMinus, yMinus, parameters); + if isempty(jumpFunction) + del_jumpy = zeros(dim_y); + del_jumpt = zeros(dim_y, 1); + else + del_jumpy = del_f_del_y(datahandle, jumpFunction, tsMinus, yPlus, parameters, h_y); + del_jumpt = del_f_del_t(datahandle, jumpFunction, tsMinus, yPlus, parameters, h_t); end - % Fix the model and set the model number for which model you want to evaluate the RHS - data.computeSensitivity.modelStage = startModel; - datahandle.setData(data); + % Evaluate the RHS at the switching point first with the model fixed on the left of the switch, then increase the model number + % and evalate the RHS with the model fixed on the left of the switch. + fminus = functionRHS(datahandle, tsMinus, yMinus, parameters); - % Calculate the remaining update matrices for the update formula until modelNum. - for i = startModel : (endModel-1) - switchingFunction = switching_functions{i}; - jumpFunction = jump_functions{i}; - ts = switches(i+1); - tsMinus = switches_left(i+1); - yMinus = y_to_switches_left(:, i+1); - yPlus = y_to_switches(:, i+1); - - h_y = fdStep_getH_y(FDstep, yMinus); - h_t = fdStep_getH_t(FDstep, tsMinus); - - % Calculate the derivatives of the switching functions w.r.t. y, t (and p if necessary) - del_sigmay = del_f_del_y(datahandle, switchingFunction, tsMinus, yMinus, parameters, h_y); - del_sigmat = del_f_del_t(datahandle, switchingFunction, tsMinus, yMinus, parameters, h_t); - diff_sigmat = del_sigmat + del_sigmay * functionRHS(datahandle, tsMinus, yMinus, parameters); - if isempty(jumpFunction) - del_jumpy = zeros(dim_y); - del_jumpt = zeros(dim_y, 1); - else - del_jumpy = del_f_del_y(datahandle, jumpFunction, tsMinus, yPlus, parameters, h_y); - del_jumpt = del_f_del_t(datahandle, jumpFunction, tsMinus, yPlus, parameters, h_t); - end + data = datahandle.getData(); + data.computeSensitivity.modelStage = data.computeSensitivity.modelStage + 1; + datahandle.setData(data); - % Evaluate the RHS at the switching point first with the model fixed on the left of the switch, then increase the model number - % and evalate the RHS with the model fixed on the left of the switch. - fminus = functionRHS(datahandle, tsMinus, yMinus, parameters); - data = datahandle.getData(); - data.computeSensitivity.modelStage = data.computeSensitivity.modelStage + 1; - datahandle.setData(data); + signatureFplus = data.SWP_detection.signature{data.computeSensitivity.modelStage}; + % Special handling for Filippov model. + if iscell(signatureFplus) + % Assume we don't start in Filippov (given since Filippov can only begin after a switch). + switchingFunction = data.SWP_detection.switchingFunction{data.computeSensitivity.modelStage - 1}; + rhsplus = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch( ... + datahandle, ... + signatureFplus{1}, ... + signatureFplus{2}, ... + switchingFunction, t, y, p); + else + rhsplus = functionRHS; + end - fplus = functionRHS(datahandle, ts, yPlus, parameters); + fplus = rhsplus(datahandle, ts, yPlus, parameters); - % Calculate the updates according to the update formula - Updates.Uy_new{i} = unit + del_jumpy + (fplus - (del_jumpt + (unit+del_jumpy)*fminus)) * del_sigmay / diff_sigmat; + % Calculate the updates according to the update formula + Updates.Uy_new{i} = unit + del_jumpy + (fplus - (del_jumpt + (unit+del_jumpy)*fminus)) * del_sigmay / diff_sigmat; - if Gp_flag - del_sigmap = del_f_del_p(datahandle, switchingFunction, tsMinus, yMinus, parameters, h_p); - if isempty(jumpFunction) - del_jumpp = zeros(dim_y, dim_p); - else - del_jumpp = del_f_del_p(datahandle, jumpFunction, tsMinus, yMinus, parameters, h_p); - end - Updates.Up_new{i} = del_jumpp + (fplus - (del_jumpt + (unit+del_jumpy)*fminus)) * del_sigmap / diff_sigmat; + if Gp_flag + del_sigmap = del_f_del_p(datahandle, switchingFunction, tsMinus, yMinus, parameters, h_p); + if isempty(jumpFunction) + del_jumpp = zeros(dim_y, dim_p); + else + del_jumpp = del_f_del_p(datahandle, jumpFunction, tsMinus, yMinus, parameters, h_p); end + Updates.Up_new{i} = del_jumpp + (fplus - (del_jumpt + (unit+del_jumpy)*fminus)) * del_sigmap / diff_sigmat; end -end \ No newline at end of file +end +end diff --git a/toolbox/internal/sensitivities/solveDisturbed_Gp.m b/toolbox/internal/sensitivities/solveDisturbed_Gp.m index 61e2da6e..60496a89 100644 --- a/toolbox/internal/sensitivities/solveDisturbed_Gp.m +++ b/toolbox/internal/sensitivities/solveDisturbed_Gp.m @@ -2,32 +2,45 @@ %SOLVEDISTURBED_GP Solve the IVP in the interval tspan with slightly disturbed parameters. % sol_original is the undisturbed solution; sols_disturbed is an array of solutions: each one has the % parameter disturbed in one component. - config = makeConfig(); - data = datahandle.getData(); +config = makeConfig(); +data = datahandle.getData(); + +signature = data.SWP_detection.signature{modelNum}; +% Special handling for Filippov model. +if iscell(signature) + % Assume we don't start in Filippov (given since Filippov can only begin after a switch). + switchingFunction = data.SWP_detection.switchingFunction{modelNum - 1}; + functionRHS_original = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch( ... + datahandle, ... + signature{1}, ... + signature{2}, ... + switchingFunction, t, y, p); +else if config.removeCtrlifForSensComputation functionRHS_original = getModelFunction(datahandle, modelNum); else functionRHS_original = data.integratorSettings.preprocessed_rhs; end - parameters = data.SWP_detection.parameters; - functionRHS_simple_END = @(t,y) functionRHS_original(datahandle, t, y, data.SWP_detection.parameters); - dim_p = data.computeSensitivity.dim_p; +end - unit_p = eye(dim_p); +parameters = data.SWP_detection.parameters; +functionRHS_simple_END = @(t,y) functionRHS_original(datahandle, t, y, data.SWP_detection.parameters); +dim_p = data.computeSensitivity.dim_p; - integrator = sensOptions.integrator; - integrator_options = sensOptions.integrator_options; +unit_p = eye(dim_p); - data.computeSensitivity.modelStage = modelNum; - datahandle.setData(data); +integrator = sensOptions.integrator; +integrator_options = sensOptions.integrator_options; - sol_original = integrator(functionRHS_simple_END, tspan, y_start, integrator_options); - sols_disturbed = cell(dim_p, 1); - for j=1:dim_p - parameters_new = parameters + h_p.*unit_p(:,j); - functionRHS_simple_disturb = @(t,y) functionRHS_original(datahandle, t, y, parameters_new); - sols_disturbed{j} = integrator(functionRHS_simple_disturb, tspan, y_start, integrator_options); - end - +data.computeSensitivity.modelStage = modelNum; +datahandle.setData(data); + +sol_original = integrator(functionRHS_simple_END, tspan, y_start, integrator_options); +sols_disturbed = cell(dim_p, 1); +for j=1:dim_p + parameters_new = parameters + h_p.*unit_p(:,j); + functionRHS_simple_disturb = @(t,y) functionRHS_original(datahandle, t, y, parameters_new); + sols_disturbed{j} = integrator(functionRHS_simple_disturb, tspan, y_start, integrator_options); end +end diff --git a/toolbox/internal/sensitivities/solveDisturbed_Gy.m b/toolbox/internal/sensitivities/solveDisturbed_Gy.m index 48fd2dbd..793c0bd3 100644 --- a/toolbox/internal/sensitivities/solveDisturbed_Gy.m +++ b/toolbox/internal/sensitivities/solveDisturbed_Gy.m @@ -2,27 +2,40 @@ %SOLVEDISTURBED_GY Solve the IVP in the interval tspan with slightly disturbed initial values. % sol_original is the undisturbed solution; sols_disturbed is an array of solutions: each one has the initial % y value disturbed in one component. - config = makeConfig(); - data = datahandle.getData(); +config = makeConfig(); +data = datahandle.getData(); + +signature = data.SWP_detection.signature{modelNum}; +% Special handling for Filippov model. +if iscell(signature) + % Assume we don't start in Filippov (given since Filippov can only begin after a switch). + switchingFunction = data.SWP_detection.switchingFunction{modelNum - 1}; + functionRHS_original = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch( ... + datahandle, ... + signature{1}, ... + signature{2}, ... + switchingFunction, t, y, p); +else if config.removeCtrlifForSensComputation functionRHS_original = getModelFunction(datahandle, modelNum); else functionRHS_original = data.integratorSettings.preprocessed_rhs; end - functionRHS_simple_END = @(t,y) functionRHS_original(datahandle, t, y, data.SWP_detection.parameters); - dim_y = data.computeSensitivity.dim_y; - unit_y = eye(dim_y); +end - integrator = sensOptions.integrator; - integrator_options = sensOptions.integrator_options; +functionRHS_simple_END = @(t,y) functionRHS_original(datahandle, t, y, data.SWP_detection.parameters); +dim_y = data.computeSensitivity.dim_y; +unit_y = eye(dim_y); - data.computeSensitivity.modelStage = modelNum; - datahandle.setData(data); +integrator = sensOptions.integrator; +integrator_options = sensOptions.integrator_options; - sol_original = integrator(functionRHS_simple_END, tspan, y_start, integrator_options); - sols_disturbed = cell(dim_y, 1); - for j=1:dim_y - sols_disturbed{j} = integrator(functionRHS_simple_END, tspan, y_start + h_y.*unit_y(:,j), integrator_options); - end -end +data.computeSensitivity.modelStage = modelNum; +datahandle.setData(data); +sol_original = integrator(functionRHS_simple_END, tspan, y_start, integrator_options); +sols_disturbed = cell(dim_y, 1); +for j=1:dim_y + sols_disturbed{j} = integrator(functionRHS_simple_END, tspan, y_start + h_y.*unit_y(:,j), integrator_options); +end +end diff --git a/toolbox/internal/sensitivities/solveVDE_Gp.m b/toolbox/internal/sensitivities/solveVDE_Gp.m index 9a656e58..5d262b73 100644 --- a/toolbox/internal/sensitivities/solveVDE_Gp.m +++ b/toolbox/internal/sensitivities/solveVDE_Gp.m @@ -1,25 +1,38 @@ function solVDE = solveVDE_Gp(datahandle, sol, tspan, modelNum, sensOptions) %SOLVEVDE_GP Solve the VDE for Gp in the interval tspan fixed to model modelNum and return the sol object - config = makeConfig(); - data = datahandle.getData(); +config = makeConfig(); +data = datahandle.getData(); - parameters = data.SWP_detection.parameters; +parameters = data.SWP_detection.parameters; + +signature = data.SWP_detection.signature{modelNum}; +% Special handling for Filippov model. +if iscell(signature) + % Assume we don't start in Filippov (given since Filippov can only begin after a switch). + switchingFunction = data.SWP_detection.switchingFunction{modelNum - 1}; + rhs = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch( ... + datahandle, ... + signature{1}, ... + signature{2}, ... + switchingFunction, t, y, p); +else if config.removeCtrlifForSensComputation rhs = getModelFunction(datahandle, modelNum); else rhs = data.integratorSettings.preprocessed_rhs; end - dim_y = data.computeSensitivity.dim_y; - dim_p = data.computeSensitivity.dim_p; +end - integrator = sensOptions.integrator; - integrator_options = sensOptions.integrator_options; +dim_y = data.computeSensitivity.dim_y; +dim_p = data.computeSensitivity.dim_p; - data.computeSensitivity.modelStage = modelNum; - datahandle.setData(data); +integrator = sensOptions.integrator; +integrator_options = sensOptions.integrator_options; - function_p_VDE = @(t,G) VDE_RHS_p(datahandle, sol, rhs, t, G, parameters, sensOptions); - initial_y = reshape(zeros(dim_y, dim_p), [], 1); - solVDE = integrator(function_p_VDE, tspan, initial_y, integrator_options); -end +data.computeSensitivity.modelStage = modelNum; +datahandle.setData(data); +function_p_VDE = @(t,G) VDE_RHS_p(datahandle, sol, rhs, t, G, parameters, sensOptions); +initial_y = reshape(zeros(dim_y, dim_p), [], 1); +solVDE = integrator(function_p_VDE, tspan, initial_y, integrator_options); +end diff --git a/toolbox/internal/sensitivities/solveVDE_Gy.m b/toolbox/internal/sensitivities/solveVDE_Gy.m index 887bf6e6..aa7f84ee 100644 --- a/toolbox/internal/sensitivities/solveVDE_Gy.m +++ b/toolbox/internal/sensitivities/solveVDE_Gy.m @@ -1,25 +1,37 @@ function solVDE = solveVDE_Gy(datahandle, sol, tspan, modelNum, sensOptions) %SOLVEVDE_GY Solve the VDE for Gy in the interval tspan fixed to model modelNum and return the sol object - config = makeConfig(); - data = datahandle.getData(); +config = makeConfig(); +data = datahandle.getData(); +signature = data.SWP_detection.signature{modelNum}; +% Special handling for Filippov model. +if iscell(signature) + % Assume we don't start in Filippov (given since Filippov can only begin after a switch). + switchingFunction = data.SWP_detection.switchingFunction{modelNum - 1}; + rhs = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch( ... + datahandle, ... + signature{1}, ... + signature{2}, ... + switchingFunction, t, y, p); +else if config.removeCtrlifForSensComputation rhs = getModelFunction(datahandle, modelNum); else rhs = data.integratorSettings.preprocessed_rhs; end - parameters = data.SWP_detection.parameters; - dim_y = data.computeSensitivity.dim_y; +end +parameters = data.SWP_detection.parameters; +dim_y = data.computeSensitivity.dim_y; - integrator = sensOptions.integrator; - integrator_options = sensOptions.integrator_options; - data.computeSensitivity.modelStage = modelNum; - datahandle.setData(data); +integrator = sensOptions.integrator; +integrator_options = sensOptions.integrator_options; - function_y_VDE = @(t,G) VDE_RHS_y(datahandle, sol, rhs, t, G, parameters, sensOptions); - initial_y = reshape(eye(dim_y), [], 1); - solVDE = integrator(function_y_VDE, tspan, initial_y, integrator_options); -end +data.computeSensitivity.modelStage = modelNum; +datahandle.setData(data); +function_y_VDE = @(t,G) VDE_RHS_y(datahandle, sol, rhs, t, G, parameters, sensOptions); +initial_y = reshape(eye(dim_y), [], 1); +solVDE = integrator(function_y_VDE, tspan, initial_y, integrator_options); +end From 296d3b2bd320a9b351ed97832146863d39fcd996 Mon Sep 17 00:00:00 2001 From: Lev Gromov Date: Mon, 17 Nov 2025 19:54:54 +0100 Subject: [PATCH 06/14] Allow bisection with same sign if close to zero via shift --- .../solving/solveODE_computeSwitchingPoint.m | 2 +- ...uteSwitchingPoint_checkSwitchingFunction.m | 46 ++++++++++++++++--- 2 files changed, 40 insertions(+), 8 deletions(-) diff --git a/toolbox/internal/solving/solveODE_computeSwitchingPoint.m b/toolbox/internal/solving/solveODE_computeSwitchingPoint.m index 60ed72fe..e9faa458 100644 --- a/toolbox/internal/solving/solveODE_computeSwitchingPoint.m +++ b/toolbox/internal/solving/solveODE_computeSwitchingPoint.m @@ -8,7 +8,7 @@ function solveODE_computeSwitchingPoint(datahandle) % preallocate for switchingIntervals -numberofswitchingfunctions = length(data.SWP_detection.switchingIndices); +numberofswitchingfunctions = length(data.SWP_detection.switchingfunctionhandles); switchingpoints = zeros(1, numberofswitchingfunctions); data.caseCtrlif = config.caseCtrlif.default; % find zero crossing of the switching function in the time interval diff --git a/toolbox/internal/solving/solveODE_computeSwitchingPoint_checkSwitchingFunction.m b/toolbox/internal/solving/solveODE_computeSwitchingPoint_checkSwitchingFunction.m index 1f077445..9e4dc11f 100644 --- a/toolbox/internal/solving/solveODE_computeSwitchingPoint_checkSwitchingFunction.m +++ b/toolbox/internal/solving/solveODE_computeSwitchingPoint_checkSwitchingFunction.m @@ -2,25 +2,57 @@ % Function that ensures that the following bisection algorithm can run % without issues. % In particular, it ensure that bisection.switchingFunction doesn't -% evaluate to zero at bisection.t1 and that it changes sign between +% evaluate to zero at bisection.t1 and that it changes sign between % bisection.t1 and bisection.t3. % % FUNCTION NEEDS REVISION. +% Sign change, everything is fine. +if sign(bisection.sw1) ~= sign(bisection.sw3) + return +end -if sign(bisection.sw1) == sign(bisection.sw3) - % no switch has occurred +tolZero = 1e-10; +% Both interval points are far from zero, so inputs are bad, nothing we can do. +if abs(bisection.sw1) >= tolZero && abs(bisection.sw3) >= tolZero % TODO: Can we be certain that the switching function doesn't have an % even number of zero crossings instead? Should probably address this % case by recursively cutting the interval and try if we can find a % subinterval [t1, t3*] with different signs on the edges. - error([ 'The switch between t=', num2str(bisection.t1), ' and t=', num2str(bisection.t3), ... - 'could not be determined. The signature did not change.\n']) + throw(noSignChangeException(bisection.t1, bisection.t3)); end +% At least one of the interval points is close to zero. +% Might be able to create a zero crossing via small shift in switching function. +shift = 2 * min(abs(bisection.sw1), abs(bisection.sw3)); +% If points positive shift down, else shift up. +if bisection.sw1 >= 0 + shift = -shift; +end +bisection.sw1 = bisection.sw1 + shift; +bisection.sw3 = bisection.sw3 + shift; +% Check sign change again after shift +if sign(bisection.sw1) == sign(bisection.sw3) + throw(noSignChangeException(bisection.t1, bisection.t3)); +end +bisection.switchingFunction = @(dh, t, x, p) bisection.switchingFunction(dh, t, x, p) + shift; +warnApplyingShift; -% Open question: What if bisection.sw1 == 0? Can we find an interval +% Open question: What if bisection.sw1 == 0? Can we find an interval % where at the edges the switching function is nonzero, % without provoking other issues? +end + +%% Exceptions +function e = noSignChangeException(t1, t2) +msg = 'Unable to determine switching point: No sign change in switching function between t=%g and t=%g.\n'; +e = MException('IFDIFF:Bisection:NoSignChange', msg, t1, t2); +end -end \ No newline at end of file +%% Warnings +function warnApplyingShift +msg = [ + 'Interval points for switching point computation have same sign, but close to zero.\n', ... + 'Applying small shift to switching function in an attempt to restore sign change.\n']; +warning('IFDIFF:Bisection:ApplyShift', msg); +end From a3a65404f0d0c54c6de6f1cb973060bcbbeceb66 Mon Sep 17 00:00:00 2001 From: Lev Gromov Date: Mon, 17 Nov 2025 19:57:42 +0100 Subject: [PATCH 07/14] Compute switching point when exiting Filippov regime --- .../solving/extendODE_filippov_regime.m | 5 +- .../internal/solving/extendODEuntilSwitch.m | 7 +- .../solving/extendODEuntilSwitch_Filippov.m | 99 +++++++++++++++++++ toolbox/solveODE.m | 47 ++++++--- 4 files changed, 137 insertions(+), 21 deletions(-) create mode 100644 toolbox/internal/solving/extendODEuntilSwitch_Filippov.m diff --git a/toolbox/internal/solving/extendODE_filippov_regime.m b/toolbox/internal/solving/extendODE_filippov_regime.m index 7bb4e3ef..2e82f52c 100644 --- a/toolbox/internal/solving/extendODE_filippov_regime.m +++ b/toolbox/internal/solving/extendODE_filippov_regime.m @@ -26,9 +26,6 @@ function extendODE_filippov_regime(datahandle) data.SWP_detection.solution_until_t3 = z; data.SWP_detection.t3 = data.SWP_detection.solution_until_t2.x(end); -% clear filippov data, serves also as indicator for inactive filippov mode -data.sliding = extendODE_filippov_regime_cleanup(data.sliding); - datahandle.setData(data); -end \ No newline at end of file +end diff --git a/toolbox/internal/solving/extendODEuntilSwitch.m b/toolbox/internal/solving/extendODEuntilSwitch.m index 28d96be0..9b626a96 100644 --- a/toolbox/internal/solving/extendODEuntilSwitch.m +++ b/toolbox/internal/solving/extendODEuntilSwitch.m @@ -51,6 +51,11 @@ function extendODEuntilSwitch(datahandle) % We can not start integrating from the old t2 because this would result in integration over a tiny % interval whose result would vanish due to limited floating point accuracy. data.SWP_detection.t2 = data.SWP_detection.t2 + baseOffset * 10^min(iter, config.switchingPointMaxPower); + if data.SWP_detection.t2 >= data.SWP_detection.t3 + % Okay, just use t3. + data.SWP_detection.t2 = data.SWP_detection.t3; + break + end datahandle.setData(data); % Check if there is a new signature. @@ -73,4 +78,4 @@ function extendODEuntilSwitch(datahandle) data.SWP_detection.switchingpoints{end + 1} = data.SWP_detection.t2; datahandle.setData(data) -end \ No newline at end of file +end diff --git a/toolbox/internal/solving/extendODEuntilSwitch_Filippov.m b/toolbox/internal/solving/extendODEuntilSwitch_Filippov.m new file mode 100644 index 00000000..abff8e75 --- /dev/null +++ b/toolbox/internal/solving/extendODEuntilSwitch_Filippov.m @@ -0,0 +1,99 @@ +function extendODEuntilSwitch_Filippov(datahandle) +%EXTENDODEUNTILSWITCH Extend ODE solution until switch and set new signature. +% EXTENDODEUNTILSWITCH(datahandle) +% +% Extend the ODE solution from the previous time point (t1) until the switching time (t2) using forced branching. +% Note that t2 has to be part of the new model (i.e. have a new signature). To achieve this, the numerically +% computed switching time (t2), which could still be part of the old model due to inaccuracies, +% may be increased by this function. +% +% INPUT +% Accepts input via datahandle providing the fields +% - SWP_detection.switch_cond_t1/t2 : switch condition before and after the numerically computed switch +% - SWP_detection.t2 : numerically computed switching time +% - SWP_detection.solution_until_t1 : ODE solution before the switch +% +% OUTPUT +% Writes results back to datahandle in the fields +% - SWP_detection.t2 : updated switching time (guaranteed to be part of a new model) +% - SWP_detection.solution_until_t2 : ODE solution until the updated switching time +% - SWP_detection.switchingpoints : new switching time appended +% - SWP_detection.signature : new signature appended + + +config = makeConfig(); +data = datahandle.getData(); + +% Check if the numerically computed switching point is part of the new model. +extendODEuntilSwitch_t1_to_t2(datahandle); +done = isFilippovExit(datahandle); + +% The numerically computed switching point may still be part of the old model due to inaccuracies. +% In that case, slightly increment the switching point and reintegrate with forced branching. +% Repeat this process until a new signature is detected. +t2FromRootFinding = data.SWP_detection.t2; +baseOffset = 16*eps(data.SWP_detection.t2); +iter = 0; +while ~done + data = datahandle.getData(); + + % Throw error if error of switching point obtained from root finding relative to tspan exceeds threshold. + % Note: No risk of division by zero since we need to make at least one integration step to detect a switch. + switchingPointError = abs((data.SWP_detection.t2 - t2FromRootFinding) ... + / (data.SWP_detection.tspan(end) - data.SWP_detection.tspan(1))); + if switchingPointError > config.switchingPointErrorThreshold + errorMsg = 'Relative error of numerically computed switching point exceeds threshold: %.16g > %.16g\n'; + error('IFDIFF:switchingPointErrorThreshold', errorMsg, switchingPointError, config.switchingPointErrorThreshold); + end + + % Increase t2 and integrate again starting from t1. + % We can not start integrating from the old t2 because this would result in integration over a tiny + % interval whose result would vanish due to limited floating point accuracy. + data.SWP_detection.t2 = data.SWP_detection.t2 + baseOffset * 10^min(iter, config.switchingPointMaxPower); + if data.SWP_detection.t2 >= data.SWP_detection.t3 + % Okay, just use t3. + data.SWP_detection.t2 = data.SWP_detection.t3; + break + end + datahandle.setData(data); + + % Check if there is a new signature. + extendODEuntilSwitch_t1_to_t2(datahandle); + done = isFilippovExit(datahandle); + + iter = iter + 1; +end + +data = datahandle.getData(); + +% Log events related to searching for the switch. +if config.debugMode + switchingPointInfo = "DEBUG: Detected switch after %d iterations at time point t=%.16g\n"; + fprintf(switchingPointInfo, iter, data.SWP_detection.t2); +end + +% Add new switching point to history. +data.SWP_detection.switchingpoints{end + 1} = data.SWP_detection.t2; + +datahandle.setData(data) +end + +%% Helpers +function updateAlpha(datahandle) +data = datahandle.getData(); +t = data.SWP_detection.t2; +x = deval(data.SWP_detection.solution_until_t2, t); +unused = data.sliding.filippov_rhs(datahandle, t, x, data.SWP_detection.parameters); %#ok +end + +function tf = isFilippovExit(datahandle) +updateAlpha(datahandle); + +data = datahandle.getData(); +[abstol, reltol] = getIntegratorTolerances(data.integratorSettings.options); +state = data.SWP_detection.solution_until_t2.x(end); +alpha_tol = min(abstol, reltol*norm(state)); +alpha = data.sliding.alpha_last; + +tf = alpha <= alpha_tol || alpha >= 1-alpha_tol; +end diff --git a/toolbox/solveODE.m b/toolbox/solveODE.m index 7bba4a6b..0f653d33 100644 --- a/toolbox/solveODE.m +++ b/toolbox/solveODE.m @@ -60,26 +60,41 @@ if checkForSwitchingIndices(datahandle) error('IFDIFF:switchInFilippov', 'Switching during or right after Filippov mode.'); end + data = datahandle.getData(); if data.SWP_detection.t2 == tspan(2) - break; % reached end of time span - else - if config.debugMode - fprintf('Left Filippov regime.\n'); - end - % Left the Filippov model and continue integration with new model: - % Need to set switching point, switching function and new signature. - % TODO: Do we need an accurate exit time of the Filippov model, e.g. for sensitivities? - % TODO: Same goes for the switching function... - % For now, just take the time where we ended up leaving and same switching function as in entry into Filippov. - data = datahandle.getData(); - data.SWP_detection.switchingpoints{end + 1} = data.SWP_detection.t2; - % TODO: This in particular seems very wrong! - data.SWP_detection.switchingFunction{end + 1} = data.SWP_detection.switchingFunction{end}; - data.SWP_detection.jumpFunction{end + 1} = []; + % clear filippov data, serves also as indicator for inactive filippov mode + data.sliding = extendODE_filippov_regime_cleanup(data.sliding); datahandle.setData(data); - solveODE_prepareNextStage(datahandle); + break; % reached end of time span + end + + % Left the Filippov model and continue integration with new model: + % Need to set switching point, switching function and new signature. + if config.debugMode + fprintf('Left Filippov regime.\n'); end + + solveODE_solution_until_t1(datahandle) + + data = datahandle.getData(); + % Switching function for entry into Filippov regime. + data.SWP_detection.switchingfunctionhandles = data.SWP_detection.switchingFunction(end); + data.SWP_detection.switchingIndices = 1; % Doesn't matter, but required for computeSwitchingPoint + % Assume no jumps in Filippov regime. + data.SWP_detection.jumpFunction{end + 1} = []; + datahandle.setData(data); + solveODE_computeSwitchingPoint(datahandle); + + % Compute actual switching point by extending solution until alpha changes. + extendODEuntilSwitch_Filippov(datahandle); + + % clear filippov data, serves also as indicator for inactive filippov mode + data = datahandle.getData(); + data.sliding = extendODE_filippov_regime_cleanup(data.sliding); + datahandle.setData(data); + + solveODE_prepareNextStage(datahandle); end % extend solution object from t2 ongoing until the next switch occurs From 4d3e57942309fa526e1c73f10895ff9c49bae369 Mon Sep 17 00:00:00 2001 From: Lev Gromov Date: Mon, 17 Nov 2025 19:59:36 +0100 Subject: [PATCH 08/14] Handle Filippov model for sensitivity update matrix --- .../internal/sensitivities/getGy_Gp_update.m | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/toolbox/internal/sensitivities/getGy_Gp_update.m b/toolbox/internal/sensitivities/getGy_Gp_update.m index 4005f5dc..bdfa7d38 100644 --- a/toolbox/internal/sensitivities/getGy_Gp_update.m +++ b/toolbox/internal/sensitivities/getGy_Gp_update.m @@ -53,10 +53,23 @@ h_y = fdStep_getH_y(FDstep, yMinus); h_t = fdStep_getH_t(FDstep, tsMinus); + signatureFminus = data.SWP_detection.signature{data.computeSensitivity.modelStage}; + % Special handling for Filippov model. + if iscell(signatureFminus) + % Assume we don't start in Filippov (given since Filippov can only begin after a switch). + switchingFunction = data.SWP_detection.switchingFunction{data.computeSensitivity.modelStage - 1}; + rhsminus = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch( ... + datahandle, ... + signatureFminus{1}, ... + signatureFminus{2}, ... + switchingFunction, t, y, p); + else + rhsminus = functionRHS; + end % Calculate the derivatives of the switching functions w.r.t. y, t (and p if necessary) del_sigmay = del_f_del_y(datahandle, switchingFunction, tsMinus, yMinus, parameters, h_y); del_sigmat = del_f_del_t(datahandle, switchingFunction, tsMinus, yMinus, parameters, h_t); - diff_sigmat = del_sigmat + del_sigmay * functionRHS(datahandle, tsMinus, yMinus, parameters); + diff_sigmat = del_sigmat + del_sigmay * rhsminus(datahandle, tsMinus, yMinus, parameters); if isempty(jumpFunction) del_jumpy = zeros(dim_y); del_jumpt = zeros(dim_y, 1); @@ -67,7 +80,8 @@ % Evaluate the RHS at the switching point first with the model fixed on the left of the switch, then increase the model number % and evalate the RHS with the model fixed on the left of the switch. - fminus = functionRHS(datahandle, tsMinus, yMinus, parameters); + + fminus = rhsminus(datahandle, tsMinus, yMinus, parameters); data = datahandle.getData(); data.computeSensitivity.modelStage = data.computeSensitivity.modelStage + 1; From 390b5a693d46485c5008637be5dce0bdae7414fc Mon Sep 17 00:00:00 2001 From: Lev Gromov Date: Mon, 17 Nov 2025 20:00:19 +0100 Subject: [PATCH 09/14] WIP: Compute sensitvities for Filippov predator-prey example --- .../examples/predatorpreyFilippov/testpprhs.m | 28 +++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/toolbox/examples/predatorpreyFilippov/testpprhs.m b/toolbox/examples/predatorpreyFilippov/testpprhs.m index 1dc27952..05739fbe 100644 --- a/toolbox/examples/predatorpreyFilippov/testpprhs.m +++ b/toolbox/examples/predatorpreyFilippov/testpprhs.m @@ -77,6 +77,24 @@ hEuler = plotit(fignum, Y_euler, 'c', namePlain(intEuler), linewidth); end +%% Sensitivity +dimy = length(x0_1); +dimp = length(p); +FDstep = generateFDstep(dimy, dimp); + +sensFunVDE = generateSensitivityFunction(datahandle, sol_ifdiff, FDstep, 'method', 'VDE', 'Gmatrices_intermediate', true); +sensFunEND = generateSensitivityFunction(datahandle, sol_ifdiff, FDstep, 'method', 'END_piecewise', 'Gmatrices_intermediate', true); + +sensVDE = sensFunVDE(X_plot); +sensEND = sensFunEND(X_plot); + +%% Plot +VDEGy11 = arrayfun(@(x) x.Gy(1,1), sensVDE); +fignum = fignum + 1; sensPlot(fignum, X_plot, VDEGy11, 'VDE'); +ENDGy11 = arrayfun(@(x) x.Gy(1,1), sensEND); +fignum = fignum + 1; sensPlot(fignum, X_plot, ENDGy11, 'END'); + + %% ANALYSIS if doErrorplot if (doIfdiff && doEuler) @@ -91,6 +109,16 @@ %% HELPERS +function sensPlot(fignum, t, G, method) + figure(fignum); clf(fignum); + semilogy(t, G, t, -G); + xlabel('t'); + ylabel('(G_y)_{11}'); + legend('G', '-G'); + title('Sens %s wrt y_1 for y_1', method); + drawnow +end + function errorPlot(fignum, x1, y1, y2, intname1, intname2) figure(fignum); clf(fignum); ydiff = calcDiff(y1, y2); From 8bc9ae60bfb420457bd28b26ee9a63170a264633 Mon Sep 17 00:00:00 2001 From: Lev Gromov Date: Sun, 30 Nov 2025 15:40:34 +0100 Subject: [PATCH 10/14] TMP(Undo before merge): Enable convexification storage, but without signatures Enabled storage of Filippov convexification data globally. Temporarily removed storage of Filippov signatures, since it does not work at the moment. Added scripts and functions for creating plots for 3d-predator-prey Filippov model for use in presentation slides. --- .../predatorpreyFilippov/exportPlotForSlide.m | 64 ++++++ .../predatorpreyFilippov/predatorPrey3D_rhs.m | 39 ++++ .../predatorPrey3D_solve.m | 208 ++++++++++++++++++ .../analyseSignature_storeSlidingInfo.m | 15 +- .../solving/cutSteps_convexification.m | 6 +- .../solving/initDatahandleFields_Filippov.m | 2 +- 6 files changed, 317 insertions(+), 17 deletions(-) create mode 100644 toolbox/examples/predatorpreyFilippov/exportPlotForSlide.m create mode 100644 toolbox/examples/predatorpreyFilippov/predatorPrey3D_rhs.m create mode 100644 toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m diff --git a/toolbox/examples/predatorpreyFilippov/exportPlotForSlide.m b/toolbox/examples/predatorpreyFilippov/exportPlotForSlide.m new file mode 100644 index 00000000..479d13db --- /dev/null +++ b/toolbox/examples/predatorpreyFilippov/exportPlotForSlide.m @@ -0,0 +1,64 @@ +function exportPlotForSlide(filepath) +set(gca, 'FontSize', 16); +set(gca, 'LineWidth', 1.2); +set(gcf, 'Color', 'w'); + +fontsize_label = 18; +removeLabels(); +%addLabelPredatorPrey(fontsize_label); +%addLabelSwitchingAlpha(fontsize_label); +addLabelSens(fontsize_label); + +%removeLegend(); +%addLegendPredatorPrey(); +addLegendSens(); + +exportgraphics(gcf, filepath, 'Resolution', 300); +end + +%% Labels +function removeLabels() +xlabel(''); +ylabel(''); +zlabel(''); +end + +function removeLabelsTwoAx() +xlabel(''); +yyaxis('left'); +ylabel(''); +yyaxis('right'); +ylabel(''); +end + +function addLabelPredatorPrey(fontsize_label) +xlabel('Räuber', 'FontSize', fontsize_label); +ylabel('Beute 2', 'FontSize', fontsize_label); +zlabel('Beute 1', 'FontSize', fontsize_label); +end + +function addLabelSwitchingAlpha(fontsize_label) +xlabel('Zeit', 'FontSize', fontsize_label) +yyaxis('left') +ylabel('Konvexitätsparameter', 'FontSize', fontsize_label) +yyaxis('right') +ylabel('Schaltfunktion', 'FontSize', fontsize_label) +end + +function addLabelSens(fontsize_label) +xlabel('Zeit', 'FontSize', fontsize_label); +ylabel('Sensitivität Räuber Anfangswert', 'FontSize', fontsize_label); +end + +%% Legend +function removeLegend() +legend('hide'); +end + +function addLegendPredatorPrey() +legend('IFDIFF', 'expl. Euler'); +end + +function addLegendSens() +legend('IFDIFF', 'expl. Euler', 'Location', 'northwest'); +end diff --git a/toolbox/examples/predatorpreyFilippov/predatorPrey3D_rhs.m b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_rhs.m new file mode 100644 index 00000000..68f3fefd --- /dev/null +++ b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_rhs.m @@ -0,0 +1,39 @@ +function dx = predatorPrey3D_rhs(~, x, p) +% 2-Prey-1-Predatory Model +% +% Source: +% Tiago Carvalho, Douglas Duarte Novaes, Luiz Fernando Goncalves: +% Sliding Shilnikov connection in Filippov-type predator-prey model +% Nonlinear Dyn (2020) 100:2973-2987 + +% Preallocate +dx = zeros(3, 1); + +% Unpack state components and parameters +p1 = x(1); % prey1 +p2 = x(2); % prey2 +P = x(3); % Predator + +r1 = p(1); +r2 = p(2); +beta1 = p(3); +beta2 = p(4); +q1 = p(5); +q2 = p(6); +m = p(7); +e = p(8); +aq = p(9); + +% Switching manifold +H = beta1*p1 - aq*beta2*p2; + +if H >= 0 + dx(1) = (r1 - beta1*P) .* p1; + dx(2) = r2*p2; + dx(3) = (e*q1*beta1*p1 - m) .* P; +else + dx(1) = r1*p1; + dx(2) = (r2 - beta2*P) .* p2; + dx(3) = (e*q2*beta2*p2 - m) .* P; +end +end diff --git a/toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m new file mode 100644 index 00000000..b9a24374 --- /dev/null +++ b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m @@ -0,0 +1,208 @@ +%% Setup +% Time span +tspan = [0, 100]; + +% Integrator +intIfdiff = @ode45; +intOptions = odeset('reltol', 1e-5, 'abstol', 1e-12, 'MaxStep', 0.5); +eulerStep = 1e-7; + +% Parameter values required for Shilnikov behavior in the paper (see RHS file). +m = 0.790; +r1 = 0.836; +e = 0.948; +q1 = 0.772; +aq = 0.660; +beta2 = 0.896; +% Some values changed for prettier plot (remains a Filippov system). +q2 = 1.500; % 1.084 (paper) +r2 = 0.300; % 0.126 (paper) +beta1 = 7.810; % 7.800 (Example plot in paper, but can be chosen arbitrarily, retaining Shilnikov behavior) +p = [r1, r2, beta1, beta2, q1, q2, m, e, aq]; + +% Initial values (used for plot in paper, but may be changed) +a = 0.286975; +x0 = [a; a; r1-r2]; + +% RHS function dx = f(t, x, p), must be implemented in separate file for IFDIFF. +rhs = @predatorPrey3D_rhs; + +% Plotting +plot_n = 10000; +plot_t = linspace(tspan(1), tspan(end), plot_n); + +nameIfdiff = sprintf('Ifdiff/%s', func2str(intIfdiff)); +nameEuler = sprintf('Euler/%.0e', eulerStep); +lwIfdiff = 3; +lwEuler = 3; +colorIfdiff = [0, 0.447, 0.741]; +colorEuler = [0.85, 0.325, 0.098]; +lsIfdiff = '-'; +lsEuler = '--'; + +% Sensitivity +wrt_y = 3; +eulerDisturbH = 1e-6; + + +%% Solve with IFDIFF +datahandle = prepareDatahandleForIntegration(rhs, 'solver', intIfdiff, 'options', intOptions); +solIfdiff = solveODE(datahandle, tspan, x0, p); + + +%% Solve with Euler +[owndir, ~] = fileparts(mfilename('fullpath')); +if contains(owndir, 'Editor_') + [owndir, ~] = fileparts(matlab.desktop.editor.getActiveFilename); +end + +euler_prefix = sprintf('sol_euler_%.0e', eulerStep); +euler_vname = 'sol_euler'; +euler_fname = fullfile(owndir, [euler_prefix, '.mat']); + +solEuler = loadEulerOrCompute(euler_fname, euler_vname, rhs, tspan, x0, p, eulerStep); + + +%% Plot solution +plot_x_ifdiff = deval(solIfdiff, plot_t); +axSol = plotSol3d([], plot_x_ifdiff, nameIfdiff, lwIfdiff, colorIfdiff, lsIfdiff); + +plot_x_euler = interp1(solEuler.x, solEuler.y', plot_t)'; +plotSol3d(axSol, plot_x_euler, nameEuler, lwEuler, colorEuler, lsEuler); + + +%% Compute Sensitivity with IFDIFF +FDstep = generateFDstep(length(x0), length(p)); +sensFunVDE = generateSensitivityFunction(datahandle, solIfdiff, FDstep, ... + 'method', 'VDE', ... + 'calcGy', true, ... + 'calcGp', false); +sensVDE = sensFunVDE(plot_t); +G = arrayfun(@(s) s.Gy(:, wrt_y), sensVDE, 'UniformOutput', false); +plot_sens_ifdiff = [G{:}]; + + +%% Compute Sensitivity with Euler +euler_prefix_disturb = sprintf('%s_disturb_wrt_y%d', euler_prefix, wrt_y); +euler_fname_disturb = fullfile(owndir, [euler_prefix_disturb, '.mat']); + +x0_disturb = x0; +x0_disturb(wrt_y) = x0_disturb(wrt_y) + eulerDisturbH; +solEulerDisturb = loadEulerOrCompute(euler_fname_disturb, euler_vname, rhs, tspan, x0_disturb, p, eulerStep); +sensEuler = (solEulerDisturb.y - solEuler.y) / eulerDisturbH; +plot_sens_euler = interp1(solEuler.x, sensEuler', plot_t)'; + + +%% Plot Sensitivity +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); +end + + +%% Plot alpha and switching function +data = datahandle.getData(); +t_alpha = data.sliding.convexification.t; +alpha = data.sliding.convexification.alpha; +ax = plotSwitchingFuncOrAlpha([], t_alpha, alpha, 'Alpha', 2, colorIfdiff, 'left', 'Alpha'); + +sw = solIfdiff.switchingFunction{1}; +sw_x = arrayfun(@(t) sw([], t, deval(solIfdiff, t), p), t_alpha); +plotSwitchingFuncOrAlpha(ax, t_alpha, sw_x, 'Switching Function', 2, colorEuler, 'right', 'Switching Function'); + +plotSwitches(ax, solIfdiff.switches) + + +%% Helpers +function sol_euler = loadEulerOrCompute(fname, vname, rhs, tspan, x0, p, step) +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('Saving result to %s for later reuse.\n', fname); +save(fname, vname); +end + +function plotSwitches(ax, switches) +xline(ax, switches, '--', 'LineWidth', 1.5, 'HandleVisibility', 'off'); +end + +function ax = plotSwitchingFuncOrAlpha(ax, t, x, name, lw, color, yyax, ylab) +if isempty(ax) + f = figure; + ax = axes(f); +end +yyaxis(ax, yyax); +hold(ax, 'on'); +plot(ax, t, x, 'DisplayName', name, 'LineWidth', lw, 'Color', color); +hold(ax, 'off'); +grid(ax, 'on'); +xlabel(ax, 't'); +ylabel(ax, ylab); +legend(ax, 'location', 'northeast'); + +yl = ylim(ax); +yl(1) = 0; +ylim(ax, yl); +end + +function ax = plotSens(ax, t, x, name, lw, color, ls) +if isempty(ax) + f = figure; + ax = axes(f); +end +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') +legend(ax, 'location', 'northeast'); +end + +function ax = plotSol3d(ax, x, name, lw, color, ls) +if isempty(ax) + f = figure; + ax = axes(f); +end +hold(ax, 'on'); +plot3(ax, x(3, :), x(2, :), x(1, :), 'DisplayName', name, 'LineWidth', lw, 'Color', color, 'LineStyle', ls); +hold(ax, 'off'); + +view(ax, [97 51]); +box(ax, 'on'); +grid(ax, 'on'); +xlabel(ax, 'Predator'); +ylabel(ax, 'Prey 2'); +zlabel(ax, 'Prey 1'); +legend(ax, 'location', 'northeast'); +end + +function sol = explEuler(rhs, tspan, x0, p, stepsize) +xdim = length(x0); +stepcount = (tspan(end)-tspan(1)) / stepsize; +sfac = 0.001; % store factor +n_out = ceil(stepcount*sfac) + 1; + +Xi = reshape(x0, [], 1); +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); + end + if ~mod(i, stepcount / 100), fprintf('.'); end +end +fprintf('\n'); +T = linspace(tspan(1), tspan(end), n_out); +sol.x = T; +sol.y = X; +end diff --git a/toolbox/internal/solving/analyseSignature_storeSlidingInfo.m b/toolbox/internal/solving/analyseSignature_storeSlidingInfo.m index 43cbe795..44b71977 100644 --- a/toolbox/internal/solving/analyseSignature_storeSlidingInfo.m +++ b/toolbox/internal/solving/analyseSignature_storeSlidingInfo.m @@ -29,18 +29,7 @@ return end -% We don't want to redundantly store the signature many times. So we store -% it once and write a pointer in form of an index into all following entries -k = length(convexification.signature_fplus); -if k==0 || ~isscalar(convexification.signature_fplus{end}) - pointer = k+1; - convexification.signature_fplus = [convexification.signature_fplus, data.sliding.signature_fplus, repmat({pointer}, 1,length(t)-1)]; - convexification.signature_fminus = [convexification.signature_fminus, data.sliding.signature_fminus, repmat({pointer}, 1,length(t)-1)]; -else - % arrays have a pointer as the last entry, so we just add copies of it - convexification.signature_fplus = [convexification.signature_fplus, repmat(convexification.signature_fplus(end), 1,length(t))]; - convexification.signature_fminus = [convexification.signature_fminus, repmat(convexification.signature_fminus(end), 1,length(t))]; -end + % t can be vector-valued, but the Filippov-RHS can only save the alpha from % the latest evaluation, so we have to recompute it instead of reading it from @@ -63,4 +52,4 @@ index = datahandle.getData().sliding.index; convexification.index = [convexification.index, repmat(index, 1,length(t))]; -end \ No newline at end of file +end diff --git a/toolbox/internal/solving/cutSteps_convexification.m b/toolbox/internal/solving/cutSteps_convexification.m index e79df427..bd39c6db 100644 --- a/toolbox/internal/solving/cutSteps_convexification.m +++ b/toolbox/internal/solving/cutSteps_convexification.m @@ -19,9 +19,9 @@ function cutSteps_convexification(datahandle, t_cut) data.sliding.convexification.t = data.sliding.convexification.t(convexification_t<=t_cut); data.sliding.convexification.index = data.sliding.convexification.index(convexification_t<=t_cut); data.sliding.convexification.alpha = data.sliding.convexification.alpha(convexification_t<=t_cut); -data.sliding.convexification.signature_fplus = data.sliding.convexification.signature_fplus(convexification_t<=t_cut); -data.sliding.convexification.signature_fminus = data.sliding.convexification.signature_fminus(convexification_t<=t_cut); +%data.sliding.convexification.signature_fplus = data.sliding.convexification.signature_fplus(convexification_t<=t_cut); +%data.sliding.convexification.signature_fminus = data.sliding.convexification.signature_fminus(convexification_t<=t_cut); datahandle.setData(data); -end \ No newline at end of file +end diff --git a/toolbox/internal/solving/initDatahandleFields_Filippov.m b/toolbox/internal/solving/initDatahandleFields_Filippov.m index 3f12be04..feda7e1d 100644 --- a/toolbox/internal/solving/initDatahandleFields_Filippov.m +++ b/toolbox/internal/solving/initDatahandleFields_Filippov.m @@ -12,7 +12,7 @@ sliding = extendODE_filippov_regime_cleanup(struct()); % if true: convexification parameters and signatures are stored for every integration step -sliding.storeSlidingInfo = false; +sliding.storeSlidingInfo = true; convexification.t = []; convexification.index = []; From b8fa5f639976d18c62dcf57f75cde55c6b03bd0e Mon Sep 17 00:00:00 2001 From: Lev Gromov Date: Tue, 27 Jan 2026 16:12:17 +0100 Subject: [PATCH 11/14] Restore old Predator-Prey example state --- .../predatorpreyFilippov/exportPlotForSlide.m | 64 ------ .../predatorpreyFilippov/predatorPrey3D_rhs.m | 39 ---- .../predatorPrey3D_solve.m | 208 ------------------ .../examples/predatorpreyFilippov/testpprhs.m | 28 --- 4 files changed, 339 deletions(-) delete mode 100644 toolbox/examples/predatorpreyFilippov/exportPlotForSlide.m delete mode 100644 toolbox/examples/predatorpreyFilippov/predatorPrey3D_rhs.m delete mode 100644 toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m diff --git a/toolbox/examples/predatorpreyFilippov/exportPlotForSlide.m b/toolbox/examples/predatorpreyFilippov/exportPlotForSlide.m deleted file mode 100644 index 479d13db..00000000 --- a/toolbox/examples/predatorpreyFilippov/exportPlotForSlide.m +++ /dev/null @@ -1,64 +0,0 @@ -function exportPlotForSlide(filepath) -set(gca, 'FontSize', 16); -set(gca, 'LineWidth', 1.2); -set(gcf, 'Color', 'w'); - -fontsize_label = 18; -removeLabels(); -%addLabelPredatorPrey(fontsize_label); -%addLabelSwitchingAlpha(fontsize_label); -addLabelSens(fontsize_label); - -%removeLegend(); -%addLegendPredatorPrey(); -addLegendSens(); - -exportgraphics(gcf, filepath, 'Resolution', 300); -end - -%% Labels -function removeLabels() -xlabel(''); -ylabel(''); -zlabel(''); -end - -function removeLabelsTwoAx() -xlabel(''); -yyaxis('left'); -ylabel(''); -yyaxis('right'); -ylabel(''); -end - -function addLabelPredatorPrey(fontsize_label) -xlabel('Räuber', 'FontSize', fontsize_label); -ylabel('Beute 2', 'FontSize', fontsize_label); -zlabel('Beute 1', 'FontSize', fontsize_label); -end - -function addLabelSwitchingAlpha(fontsize_label) -xlabel('Zeit', 'FontSize', fontsize_label) -yyaxis('left') -ylabel('Konvexitätsparameter', 'FontSize', fontsize_label) -yyaxis('right') -ylabel('Schaltfunktion', 'FontSize', fontsize_label) -end - -function addLabelSens(fontsize_label) -xlabel('Zeit', 'FontSize', fontsize_label); -ylabel('Sensitivität Räuber Anfangswert', 'FontSize', fontsize_label); -end - -%% Legend -function removeLegend() -legend('hide'); -end - -function addLegendPredatorPrey() -legend('IFDIFF', 'expl. Euler'); -end - -function addLegendSens() -legend('IFDIFF', 'expl. Euler', 'Location', 'northwest'); -end diff --git a/toolbox/examples/predatorpreyFilippov/predatorPrey3D_rhs.m b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_rhs.m deleted file mode 100644 index 68f3fefd..00000000 --- a/toolbox/examples/predatorpreyFilippov/predatorPrey3D_rhs.m +++ /dev/null @@ -1,39 +0,0 @@ -function dx = predatorPrey3D_rhs(~, x, p) -% 2-Prey-1-Predatory Model -% -% Source: -% Tiago Carvalho, Douglas Duarte Novaes, Luiz Fernando Goncalves: -% Sliding Shilnikov connection in Filippov-type predator-prey model -% Nonlinear Dyn (2020) 100:2973-2987 - -% Preallocate -dx = zeros(3, 1); - -% Unpack state components and parameters -p1 = x(1); % prey1 -p2 = x(2); % prey2 -P = x(3); % Predator - -r1 = p(1); -r2 = p(2); -beta1 = p(3); -beta2 = p(4); -q1 = p(5); -q2 = p(6); -m = p(7); -e = p(8); -aq = p(9); - -% Switching manifold -H = beta1*p1 - aq*beta2*p2; - -if H >= 0 - dx(1) = (r1 - beta1*P) .* p1; - dx(2) = r2*p2; - dx(3) = (e*q1*beta1*p1 - m) .* P; -else - dx(1) = r1*p1; - dx(2) = (r2 - beta2*P) .* p2; - dx(3) = (e*q2*beta2*p2 - m) .* P; -end -end diff --git a/toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m deleted file mode 100644 index b9a24374..00000000 --- a/toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m +++ /dev/null @@ -1,208 +0,0 @@ -%% Setup -% Time span -tspan = [0, 100]; - -% Integrator -intIfdiff = @ode45; -intOptions = odeset('reltol', 1e-5, 'abstol', 1e-12, 'MaxStep', 0.5); -eulerStep = 1e-7; - -% Parameter values required for Shilnikov behavior in the paper (see RHS file). -m = 0.790; -r1 = 0.836; -e = 0.948; -q1 = 0.772; -aq = 0.660; -beta2 = 0.896; -% Some values changed for prettier plot (remains a Filippov system). -q2 = 1.500; % 1.084 (paper) -r2 = 0.300; % 0.126 (paper) -beta1 = 7.810; % 7.800 (Example plot in paper, but can be chosen arbitrarily, retaining Shilnikov behavior) -p = [r1, r2, beta1, beta2, q1, q2, m, e, aq]; - -% Initial values (used for plot in paper, but may be changed) -a = 0.286975; -x0 = [a; a; r1-r2]; - -% RHS function dx = f(t, x, p), must be implemented in separate file for IFDIFF. -rhs = @predatorPrey3D_rhs; - -% Plotting -plot_n = 10000; -plot_t = linspace(tspan(1), tspan(end), plot_n); - -nameIfdiff = sprintf('Ifdiff/%s', func2str(intIfdiff)); -nameEuler = sprintf('Euler/%.0e', eulerStep); -lwIfdiff = 3; -lwEuler = 3; -colorIfdiff = [0, 0.447, 0.741]; -colorEuler = [0.85, 0.325, 0.098]; -lsIfdiff = '-'; -lsEuler = '--'; - -% Sensitivity -wrt_y = 3; -eulerDisturbH = 1e-6; - - -%% Solve with IFDIFF -datahandle = prepareDatahandleForIntegration(rhs, 'solver', intIfdiff, 'options', intOptions); -solIfdiff = solveODE(datahandle, tspan, x0, p); - - -%% Solve with Euler -[owndir, ~] = fileparts(mfilename('fullpath')); -if contains(owndir, 'Editor_') - [owndir, ~] = fileparts(matlab.desktop.editor.getActiveFilename); -end - -euler_prefix = sprintf('sol_euler_%.0e', eulerStep); -euler_vname = 'sol_euler'; -euler_fname = fullfile(owndir, [euler_prefix, '.mat']); - -solEuler = loadEulerOrCompute(euler_fname, euler_vname, rhs, tspan, x0, p, eulerStep); - - -%% Plot solution -plot_x_ifdiff = deval(solIfdiff, plot_t); -axSol = plotSol3d([], plot_x_ifdiff, nameIfdiff, lwIfdiff, colorIfdiff, lsIfdiff); - -plot_x_euler = interp1(solEuler.x, solEuler.y', plot_t)'; -plotSol3d(axSol, plot_x_euler, nameEuler, lwEuler, colorEuler, lsEuler); - - -%% Compute Sensitivity with IFDIFF -FDstep = generateFDstep(length(x0), length(p)); -sensFunVDE = generateSensitivityFunction(datahandle, solIfdiff, FDstep, ... - 'method', 'VDE', ... - 'calcGy', true, ... - 'calcGp', false); -sensVDE = sensFunVDE(plot_t); -G = arrayfun(@(s) s.Gy(:, wrt_y), sensVDE, 'UniformOutput', false); -plot_sens_ifdiff = [G{:}]; - - -%% Compute Sensitivity with Euler -euler_prefix_disturb = sprintf('%s_disturb_wrt_y%d', euler_prefix, wrt_y); -euler_fname_disturb = fullfile(owndir, [euler_prefix_disturb, '.mat']); - -x0_disturb = x0; -x0_disturb(wrt_y) = x0_disturb(wrt_y) + eulerDisturbH; -solEulerDisturb = loadEulerOrCompute(euler_fname_disturb, euler_vname, rhs, tspan, x0_disturb, p, eulerStep); -sensEuler = (solEulerDisturb.y - solEuler.y) / eulerDisturbH; -plot_sens_euler = interp1(solEuler.x, sensEuler', plot_t)'; - - -%% Plot Sensitivity -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); -end - - -%% Plot alpha and switching function -data = datahandle.getData(); -t_alpha = data.sliding.convexification.t; -alpha = data.sliding.convexification.alpha; -ax = plotSwitchingFuncOrAlpha([], t_alpha, alpha, 'Alpha', 2, colorIfdiff, 'left', 'Alpha'); - -sw = solIfdiff.switchingFunction{1}; -sw_x = arrayfun(@(t) sw([], t, deval(solIfdiff, t), p), t_alpha); -plotSwitchingFuncOrAlpha(ax, t_alpha, sw_x, 'Switching Function', 2, colorEuler, 'right', 'Switching Function'); - -plotSwitches(ax, solIfdiff.switches) - - -%% Helpers -function sol_euler = loadEulerOrCompute(fname, vname, rhs, tspan, x0, p, step) -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('Saving result to %s for later reuse.\n', fname); -save(fname, vname); -end - -function plotSwitches(ax, switches) -xline(ax, switches, '--', 'LineWidth', 1.5, 'HandleVisibility', 'off'); -end - -function ax = plotSwitchingFuncOrAlpha(ax, t, x, name, lw, color, yyax, ylab) -if isempty(ax) - f = figure; - ax = axes(f); -end -yyaxis(ax, yyax); -hold(ax, 'on'); -plot(ax, t, x, 'DisplayName', name, 'LineWidth', lw, 'Color', color); -hold(ax, 'off'); -grid(ax, 'on'); -xlabel(ax, 't'); -ylabel(ax, ylab); -legend(ax, 'location', 'northeast'); - -yl = ylim(ax); -yl(1) = 0; -ylim(ax, yl); -end - -function ax = plotSens(ax, t, x, name, lw, color, ls) -if isempty(ax) - f = figure; - ax = axes(f); -end -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') -legend(ax, 'location', 'northeast'); -end - -function ax = plotSol3d(ax, x, name, lw, color, ls) -if isempty(ax) - f = figure; - ax = axes(f); -end -hold(ax, 'on'); -plot3(ax, x(3, :), x(2, :), x(1, :), 'DisplayName', name, 'LineWidth', lw, 'Color', color, 'LineStyle', ls); -hold(ax, 'off'); - -view(ax, [97 51]); -box(ax, 'on'); -grid(ax, 'on'); -xlabel(ax, 'Predator'); -ylabel(ax, 'Prey 2'); -zlabel(ax, 'Prey 1'); -legend(ax, 'location', 'northeast'); -end - -function sol = explEuler(rhs, tspan, x0, p, stepsize) -xdim = length(x0); -stepcount = (tspan(end)-tspan(1)) / stepsize; -sfac = 0.001; % store factor -n_out = ceil(stepcount*sfac) + 1; - -Xi = reshape(x0, [], 1); -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); - end - if ~mod(i, stepcount / 100), fprintf('.'); end -end -fprintf('\n'); -T = linspace(tspan(1), tspan(end), n_out); -sol.x = T; -sol.y = X; -end diff --git a/toolbox/examples/predatorpreyFilippov/testpprhs.m b/toolbox/examples/predatorpreyFilippov/testpprhs.m index 05739fbe..1dc27952 100644 --- a/toolbox/examples/predatorpreyFilippov/testpprhs.m +++ b/toolbox/examples/predatorpreyFilippov/testpprhs.m @@ -77,24 +77,6 @@ hEuler = plotit(fignum, Y_euler, 'c', namePlain(intEuler), linewidth); end -%% Sensitivity -dimy = length(x0_1); -dimp = length(p); -FDstep = generateFDstep(dimy, dimp); - -sensFunVDE = generateSensitivityFunction(datahandle, sol_ifdiff, FDstep, 'method', 'VDE', 'Gmatrices_intermediate', true); -sensFunEND = generateSensitivityFunction(datahandle, sol_ifdiff, FDstep, 'method', 'END_piecewise', 'Gmatrices_intermediate', true); - -sensVDE = sensFunVDE(X_plot); -sensEND = sensFunEND(X_plot); - -%% Plot -VDEGy11 = arrayfun(@(x) x.Gy(1,1), sensVDE); -fignum = fignum + 1; sensPlot(fignum, X_plot, VDEGy11, 'VDE'); -ENDGy11 = arrayfun(@(x) x.Gy(1,1), sensEND); -fignum = fignum + 1; sensPlot(fignum, X_plot, ENDGy11, 'END'); - - %% ANALYSIS if doErrorplot if (doIfdiff && doEuler) @@ -109,16 +91,6 @@ %% HELPERS -function sensPlot(fignum, t, G, method) - figure(fignum); clf(fignum); - semilogy(t, G, t, -G); - xlabel('t'); - ylabel('(G_y)_{11}'); - legend('G', '-G'); - title('Sens %s wrt y_1 for y_1', method); - drawnow -end - function errorPlot(fignum, x1, y1, y2, intname1, intname2) figure(fignum); clf(fignum); ydiff = calcDiff(y1, y2); From 45f5daf576c669111a335b1d4a256f1ca7c85d51 Mon Sep 17 00:00:00 2001 From: Lev Gromov Date: Tue, 27 Jan 2026 17:22:07 +0100 Subject: [PATCH 12/14] Simplify retrieval of submodel RHS in sensitivity computation --- .../internal/sensitivities/getGy_Gp_update.m | 41 +++------------- .../internal/sensitivities/getModelFunction.m | 26 ---------- .../sensitivities/getRhsFromModelNum.m | 48 +++++++++++++++++++ .../sensitivities/solveDisturbed_Gp.m | 19 +------- .../sensitivities/solveDisturbed_Gy.m | 19 +------- toolbox/internal/sensitivities/solveVDE_Gp.m | 19 +------- toolbox/internal/sensitivities/solveVDE_Gy.m | 19 +------- 7 files changed, 59 insertions(+), 132 deletions(-) delete mode 100644 toolbox/internal/sensitivities/getModelFunction.m create mode 100644 toolbox/internal/sensitivities/getRhsFromModelNum.m diff --git a/toolbox/internal/sensitivities/getGy_Gp_update.m b/toolbox/internal/sensitivities/getGy_Gp_update.m index bdfa7d38..c053ff5b 100644 --- a/toolbox/internal/sensitivities/getGy_Gp_update.m +++ b/toolbox/internal/sensitivities/getGy_Gp_update.m @@ -21,7 +21,6 @@ y_to_switches_left = data.computeSensitivity.y_to_switches_left; switching_functions = data.SWP_detection.switchingFunction; jump_functions = data.SWP_detection.jumpFunction; -functionRHS = data.integratorSettings.preprocessed_rhs; FDstep = options.FDstep; @@ -41,6 +40,7 @@ data.computeSensitivity.modelStage = startModel; datahandle.setData(data); +rhs = getRhsFromModelNum(datahandle, startModel); % Calculate the remaining update matrices for the update formula until modelNum. for i = startModel : (endModel-1) switchingFunction = switching_functions{i}; @@ -53,23 +53,10 @@ h_y = fdStep_getH_y(FDstep, yMinus); h_t = fdStep_getH_t(FDstep, tsMinus); - signatureFminus = data.SWP_detection.signature{data.computeSensitivity.modelStage}; - % Special handling for Filippov model. - if iscell(signatureFminus) - % Assume we don't start in Filippov (given since Filippov can only begin after a switch). - switchingFunction = data.SWP_detection.switchingFunction{data.computeSensitivity.modelStage - 1}; - rhsminus = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch( ... - datahandle, ... - signatureFminus{1}, ... - signatureFminus{2}, ... - switchingFunction, t, y, p); - else - rhsminus = functionRHS; - end % Calculate the derivatives of the switching functions w.r.t. y, t (and p if necessary) del_sigmay = del_f_del_y(datahandle, switchingFunction, tsMinus, yMinus, parameters, h_y); del_sigmat = del_f_del_t(datahandle, switchingFunction, tsMinus, yMinus, parameters, h_t); - diff_sigmat = del_sigmat + del_sigmay * rhsminus(datahandle, tsMinus, yMinus, parameters); + diff_sigmat = del_sigmat + del_sigmay * rhs(datahandle, tsMinus, yMinus, parameters); if isempty(jumpFunction) del_jumpy = zeros(dim_y); del_jumpt = zeros(dim_y, 1); @@ -80,29 +67,15 @@ % Evaluate the RHS at the switching point first with the model fixed on the left of the switch, then increase the model number % and evalate the RHS with the model fixed on the left of the switch. + fminus = rhs(datahandle, tsMinus, yMinus, parameters); - fminus = rhsminus(datahandle, tsMinus, yMinus, parameters); - + % If the RHS has ctrlifs, this ensures that the correct branching is applied. data = datahandle.getData(); data.computeSensitivity.modelStage = data.computeSensitivity.modelStage + 1; datahandle.setData(data); - - - signatureFplus = data.SWP_detection.signature{data.computeSensitivity.modelStage}; - % Special handling for Filippov model. - if iscell(signatureFplus) - % Assume we don't start in Filippov (given since Filippov can only begin after a switch). - switchingFunction = data.SWP_detection.switchingFunction{data.computeSensitivity.modelStage - 1}; - rhsplus = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch( ... - datahandle, ... - signatureFplus{1}, ... - signatureFplus{2}, ... - switchingFunction, t, y, p); - else - rhsplus = functionRHS; - end - - fplus = rhsplus(datahandle, ts, yPlus, parameters); + % Get RHS for next submodel (will correspond to fminus in the next iteration) + rhs = getRhsFromModelNum(datahandle, i+1); + fplus = rhs(datahandle, ts, yPlus, parameters); % Calculate the updates according to the update formula Updates.Uy_new{i} = unit + del_jumpy + (fplus - (del_jumpt + (unit+del_jumpy)*fminus)) * del_sigmay / diff_sigmat; diff --git a/toolbox/internal/sensitivities/getModelFunction.m b/toolbox/internal/sensitivities/getModelFunction.m deleted file mode 100644 index 44d803b5..00000000 --- a/toolbox/internal/sensitivities/getModelFunction.m +++ /dev/null @@ -1,26 +0,0 @@ -function functionHandle = getModelFunction(datahandle, modelStage) -%functionHandle = GETMODELFUNCTION(datahandle, modelStage) -% -%Find an existing model function or create a new one based on a model stage. -% -%INPUT: -% datahandle - Contains signature corresponding to the model stage. -% struct -% -% modelStage - Number of the model as encountered during initial computation of solution. -% positive integer -% -%OUTPUT: -% functionHandle - Handle to the requested model function. -% function_handle - -data = datahandle.getData(); -factory = data.codeGen.modelFactory; - -signature = data.SWP_detection.signature{modelStage}; - -[functionHandle, collisionIndex] = factory.findExisting(signature); -if isempty(functionHandle) - functionHandle = factory.createNew(signature, collisionIndex); -end -end diff --git a/toolbox/internal/sensitivities/getRhsFromModelNum.m b/toolbox/internal/sensitivities/getRhsFromModelNum.m new file mode 100644 index 00000000..eda9de82 --- /dev/null +++ b/toolbox/internal/sensitivities/getRhsFromModelNum.m @@ -0,0 +1,48 @@ +function rhs = getRhsFromModelNum(datahandle, modelNum) +%rhs = GETRHSFROMMODELNUM(datahandle, modelNum) +% +%Return handle to RHS for a particular (potentially Filippov) submodel. +% +%If submodel RHS functions are configured to be exported, this function will attempt +%to find an existing model function or create a new one based on the submodel signature. +% +%INPUT: +% datahandle - Contains stored information about RHS. +% struct +% +% modelNum - Number of the submodel in order of occurrence during initial solution via solveODE. +% positive integer +% +%OUTPUT: +% rhs - Handle to the submodel function corresponding to the signature. +% function_handle + +config = makeConfig(); +data = datahandle.getData(); + +signature = data.SWP_detection.signature{modelNum}; +% Check for Filippov case (multiple signatures) +if iscell(signature) + % TODO: Add option to also export Filippov RHS source code without ctrlifs. + % Assume we don't start in Filippov (guaranteed since Filippov can only begin after a switch). + switchingFunction = data.SWP_detection.switchingFunction{modelNum - 1}; + rhs = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch( ... + datahandle, ... + signature{1}, ... + signature{2}, ... + switchingFunction, t, y, p); + return +end + +% Get transformed RHS source code without Ctrlifs +if config.removeCtrlifForSensComputation + factory = data.codeGen.modelFactory; + [rhs, collisionIndex] = factory.findExisting(signature); + if isempty(rhs) + rhs = factory.createNew(signature, collisionIndex); + end + return +end + +% Default case, just return the preprocessed RHS +rhs = data.integratorSettings.preprocessed_rhs; diff --git a/toolbox/internal/sensitivities/solveDisturbed_Gp.m b/toolbox/internal/sensitivities/solveDisturbed_Gp.m index 60496a89..b3aa3857 100644 --- a/toolbox/internal/sensitivities/solveDisturbed_Gp.m +++ b/toolbox/internal/sensitivities/solveDisturbed_Gp.m @@ -2,26 +2,9 @@ %SOLVEDISTURBED_GP Solve the IVP in the interval tspan with slightly disturbed parameters. % sol_original is the undisturbed solution; sols_disturbed is an array of solutions: each one has the % parameter disturbed in one component. -config = makeConfig(); data = datahandle.getData(); -signature = data.SWP_detection.signature{modelNum}; -% Special handling for Filippov model. -if iscell(signature) - % Assume we don't start in Filippov (given since Filippov can only begin after a switch). - switchingFunction = data.SWP_detection.switchingFunction{modelNum - 1}; - functionRHS_original = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch( ... - datahandle, ... - signature{1}, ... - signature{2}, ... - switchingFunction, t, y, p); -else - if config.removeCtrlifForSensComputation - functionRHS_original = getModelFunction(datahandle, modelNum); - else - functionRHS_original = data.integratorSettings.preprocessed_rhs; - end -end +functionRHS_original = getRhsFromModelNum(datahandle, modelNum); parameters = data.SWP_detection.parameters; functionRHS_simple_END = @(t,y) functionRHS_original(datahandle, t, y, data.SWP_detection.parameters); diff --git a/toolbox/internal/sensitivities/solveDisturbed_Gy.m b/toolbox/internal/sensitivities/solveDisturbed_Gy.m index 793c0bd3..1a04855b 100644 --- a/toolbox/internal/sensitivities/solveDisturbed_Gy.m +++ b/toolbox/internal/sensitivities/solveDisturbed_Gy.m @@ -2,26 +2,9 @@ %SOLVEDISTURBED_GY Solve the IVP in the interval tspan with slightly disturbed initial values. % sol_original is the undisturbed solution; sols_disturbed is an array of solutions: each one has the initial % y value disturbed in one component. -config = makeConfig(); data = datahandle.getData(); -signature = data.SWP_detection.signature{modelNum}; -% Special handling for Filippov model. -if iscell(signature) - % Assume we don't start in Filippov (given since Filippov can only begin after a switch). - switchingFunction = data.SWP_detection.switchingFunction{modelNum - 1}; - functionRHS_original = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch( ... - datahandle, ... - signature{1}, ... - signature{2}, ... - switchingFunction, t, y, p); -else - if config.removeCtrlifForSensComputation - functionRHS_original = getModelFunction(datahandle, modelNum); - else - functionRHS_original = data.integratorSettings.preprocessed_rhs; - end -end +functionRHS_original = getRhsFromModelNum(datahandle, modelNum); functionRHS_simple_END = @(t,y) functionRHS_original(datahandle, t, y, data.SWP_detection.parameters); dim_y = data.computeSensitivity.dim_y; diff --git a/toolbox/internal/sensitivities/solveVDE_Gp.m b/toolbox/internal/sensitivities/solveVDE_Gp.m index 5d262b73..1da1239b 100644 --- a/toolbox/internal/sensitivities/solveVDE_Gp.m +++ b/toolbox/internal/sensitivities/solveVDE_Gp.m @@ -1,27 +1,10 @@ function solVDE = solveVDE_Gp(datahandle, sol, tspan, modelNum, sensOptions) %SOLVEVDE_GP Solve the VDE for Gp in the interval tspan fixed to model modelNum and return the sol object -config = makeConfig(); data = datahandle.getData(); parameters = data.SWP_detection.parameters; -signature = data.SWP_detection.signature{modelNum}; -% Special handling for Filippov model. -if iscell(signature) - % Assume we don't start in Filippov (given since Filippov can only begin after a switch). - switchingFunction = data.SWP_detection.switchingFunction{modelNum - 1}; - rhs = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch( ... - datahandle, ... - signature{1}, ... - signature{2}, ... - switchingFunction, t, y, p); -else - if config.removeCtrlifForSensComputation - rhs = getModelFunction(datahandle, modelNum); - else - rhs = data.integratorSettings.preprocessed_rhs; - end -end +rhs = getRhsFromModelNum(datahandle, modelNum); dim_y = data.computeSensitivity.dim_y; dim_p = data.computeSensitivity.dim_p; diff --git a/toolbox/internal/sensitivities/solveVDE_Gy.m b/toolbox/internal/sensitivities/solveVDE_Gy.m index aa7f84ee..b297acc9 100644 --- a/toolbox/internal/sensitivities/solveVDE_Gy.m +++ b/toolbox/internal/sensitivities/solveVDE_Gy.m @@ -1,25 +1,8 @@ function solVDE = solveVDE_Gy(datahandle, sol, tspan, modelNum, sensOptions) %SOLVEVDE_GY Solve the VDE for Gy in the interval tspan fixed to model modelNum and return the sol object -config = makeConfig(); data = datahandle.getData(); -signature = data.SWP_detection.signature{modelNum}; -% Special handling for Filippov model. -if iscell(signature) - % Assume we don't start in Filippov (given since Filippov can only begin after a switch). - switchingFunction = data.SWP_detection.switchingFunction{modelNum - 1}; - rhs = @(datahandle, t, y, p) slidingFilippovRHS_oneSwitch( ... - datahandle, ... - signature{1}, ... - signature{2}, ... - switchingFunction, t, y, p); -else - if config.removeCtrlifForSensComputation - rhs = getModelFunction(datahandle, modelNum); - else - rhs = data.integratorSettings.preprocessed_rhs; - end -end +rhs = getRhsFromModelNum(datahandle, modelNum); parameters = data.SWP_detection.parameters; dim_y = data.computeSensitivity.dim_y; From 1c6690a15bfb38b8fc65d13583acd15ecf53ab79 Mon Sep 17 00:00:00 2001 From: Lev Gromov Date: Tue, 27 Jan 2026 17:57:46 +0100 Subject: [PATCH 13/14] Remove redundant storage of signature in convexification --- .../internal/solving/analyseSignature_storeSlidingInfo.m | 2 -- toolbox/internal/solving/cutSteps_convexification.m | 9 ++++----- toolbox/internal/solving/initDatahandleFields_Filippov.m | 4 +--- 3 files changed, 5 insertions(+), 10 deletions(-) diff --git a/toolbox/internal/solving/analyseSignature_storeSlidingInfo.m b/toolbox/internal/solving/analyseSignature_storeSlidingInfo.m index 44b71977..7c301d6b 100644 --- a/toolbox/internal/solving/analyseSignature_storeSlidingInfo.m +++ b/toolbox/internal/solving/analyseSignature_storeSlidingInfo.m @@ -24,8 +24,6 @@ % Not in filippov mode, add NaN/entry entries. convexification.index = [convexification.index, NaN(1,length(t))]; convexification.alpha = [convexification.alpha, NaN(1,length(t))]; - convexification.signature_fplus = [convexification.signature_fplus, cell(1,length(t))]; % write empty cells - convexification.signature_fminus = [convexification.signature_fminus, cell(1,length(t))]; return end diff --git a/toolbox/internal/solving/cutSteps_convexification.m b/toolbox/internal/solving/cutSteps_convexification.m index bd39c6db..d3f50666 100644 --- a/toolbox/internal/solving/cutSteps_convexification.m +++ b/toolbox/internal/solving/cutSteps_convexification.m @@ -16,11 +16,10 @@ function cutSteps_convexification(datahandle, t_cut) data = datahandle.getData(); convexification_t = data.sliding.convexification.t; -data.sliding.convexification.t = data.sliding.convexification.t(convexification_t<=t_cut); -data.sliding.convexification.index = data.sliding.convexification.index(convexification_t<=t_cut); -data.sliding.convexification.alpha = data.sliding.convexification.alpha(convexification_t<=t_cut); -%data.sliding.convexification.signature_fplus = data.sliding.convexification.signature_fplus(convexification_t<=t_cut); -%data.sliding.convexification.signature_fminus = data.sliding.convexification.signature_fminus(convexification_t<=t_cut); +mask = convexification_t <= t_cut; +data.sliding.convexification.t = data.sliding.convexification.t(mask); +data.sliding.convexification.index = data.sliding.convexification.index(mask); +data.sliding.convexification.alpha = data.sliding.convexification.alpha(mask); datahandle.setData(data); diff --git a/toolbox/internal/solving/initDatahandleFields_Filippov.m b/toolbox/internal/solving/initDatahandleFields_Filippov.m index feda7e1d..4972f940 100644 --- a/toolbox/internal/solving/initDatahandleFields_Filippov.m +++ b/toolbox/internal/solving/initDatahandleFields_Filippov.m @@ -12,12 +12,10 @@ sliding = extendODE_filippov_regime_cleanup(struct()); % if true: convexification parameters and signatures are stored for every integration step -sliding.storeSlidingInfo = true; +sliding.storeSlidingInfo = false; convexification.t = []; convexification.index = []; convexification.alpha = []; -convexification.signature_fminus = {}; -convexification.signature_fplus = {}; sliding.convexification = convexification; end From 67c003bdcc94556bdcf75d44556e3aae3a5c32fa Mon Sep 17 00:00:00 2001 From: Lev Gromov Date: Tue, 27 Jan 2026 18:56:44 +0100 Subject: [PATCH 14/14] Extend predator-prey example with alpha tracking and sensitivities --- toolbox/examples/predatorpreyFilippov/pprhs.m | 44 ---- .../predatorpreyFilippov/predatorPrey3D_rhs.m | 39 ++++ .../predatorPrey3D_solve.m | 213 ++++++++++++++++++ .../examples/predatorpreyFilippov/testpprhs.m | 162 ------------- .../rhsPreprocessing/createDatahandle.m | 2 + .../internal/solving/initDatahandleFields.m | 2 +- .../solving/initDatahandleFields_Filippov.m | 9 +- .../solving/solveODE_setFilippovRHS.m | 2 - toolbox/makeConfig.m | 3 + 9 files changed, 261 insertions(+), 215 deletions(-) delete mode 100644 toolbox/examples/predatorpreyFilippov/pprhs.m create mode 100644 toolbox/examples/predatorpreyFilippov/predatorPrey3D_rhs.m create mode 100644 toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m delete mode 100644 toolbox/examples/predatorpreyFilippov/testpprhs.m diff --git a/toolbox/examples/predatorpreyFilippov/pprhs.m b/toolbox/examples/predatorpreyFilippov/pprhs.m deleted file mode 100644 index e1b1a8ca..00000000 --- a/toolbox/examples/predatorpreyFilippov/pprhs.m +++ /dev/null @@ -1,44 +0,0 @@ -function dx = pprhs(t,x,p) -% Predator Prey Model -% -% Source: -% Tiago Carvalho, Douglas Duarte Novaes, Luiz Fernando Goncalves: -% Sliding Shilnikov connection in Filippov-type predator-prey model -% Nonlinear Dyn (2020) 100:2973-2987 -% - -% preallocate output -dx = [0 ; 0 ; 0]; - -% x = (p1, p2, P) % prey, prey2, predator -p1 = x(1); -p2 = x(2); -P = x(3); - -% p = (r1, r2, beta1, beta2, q1, q2, m, e, aq) % parameters -r1 = p(1); -r2 = p(2); -beta1 = p(3); -beta2 = p(4); -q1 = p(5); -q2 = p(6); -m = p(7); -e = p(8); -aq = p(9); - -% switching manifold -H = beta1*p1 - aq*beta2*p2; - -% switched rhs -if H >= 0 - dx(1) = (r1 - beta1*P) .* p1; - dx(2) = r2*p2; - dx(3) = (e*q1*beta1*p1 - m) .* P; -else - dx(1) = r1*p1; - dx(2) = (r2 - beta2*P) .* p2; - dx(3) = (e*q2*beta2*p2 - m) .* P; -end - - -end \ No newline at end of file diff --git a/toolbox/examples/predatorpreyFilippov/predatorPrey3D_rhs.m b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_rhs.m new file mode 100644 index 00000000..68f3fefd --- /dev/null +++ b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_rhs.m @@ -0,0 +1,39 @@ +function dx = predatorPrey3D_rhs(~, x, p) +% 2-Prey-1-Predatory Model +% +% Source: +% Tiago Carvalho, Douglas Duarte Novaes, Luiz Fernando Goncalves: +% Sliding Shilnikov connection in Filippov-type predator-prey model +% Nonlinear Dyn (2020) 100:2973-2987 + +% Preallocate +dx = zeros(3, 1); + +% Unpack state components and parameters +p1 = x(1); % prey1 +p2 = x(2); % prey2 +P = x(3); % Predator + +r1 = p(1); +r2 = p(2); +beta1 = p(3); +beta2 = p(4); +q1 = p(5); +q2 = p(6); +m = p(7); +e = p(8); +aq = p(9); + +% Switching manifold +H = beta1*p1 - aq*beta2*p2; + +if H >= 0 + dx(1) = (r1 - beta1*P) .* p1; + dx(2) = r2*p2; + dx(3) = (e*q1*beta1*p1 - m) .* P; +else + dx(1) = r1*p1; + dx(2) = (r2 - beta2*P) .* p2; + dx(3) = (e*q2*beta2*p2 - m) .* P; +end +end diff --git a/toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m new file mode 100644 index 00000000..0531bbdb --- /dev/null +++ b/toolbox/examples/predatorpreyFilippov/predatorPrey3D_solve.m @@ -0,0 +1,213 @@ +%% Setup +% Time span +tspan = [0, 100]; + +% Integrator +intIfdiff = @ode45; +intOptions = odeset('reltol', 1e-5, 'abstol', 1e-12, 'MaxStep', 0.5); +eulerStep = 1e-7; + +% Parameter values required for Shilnikov behavior in the paper (see RHS file). +m = 0.790; +r1 = 0.836; +e = 0.948; +q1 = 0.772; +aq = 0.660; +beta2 = 0.896; +% Some values changed for prettier plot (remains a Filippov system). +q2 = 1.500; % 1.084 (paper) +r2 = 0.300; % 0.126 (paper) +beta1 = 7.810; % 7.800 (Example plot in paper, but can be chosen arbitrarily, retaining Shilnikov behavior) +p = [r1, r2, beta1, beta2, q1, q2, m, e, aq]; + +% Initial values (used for plot in paper, but may be changed) +a = 0.286975; +x0 = [a; a; r1-r2]; + +% RHS function dx = f(t, x, p), must be implemented in separate file for IFDIFF. +rhs = @predatorPrey3D_rhs; + +% Plotting +plot_n = 10000; +plot_t = linspace(tspan(1), tspan(end), plot_n); + +nameIfdiff = sprintf('Ifdiff/%s', func2str(intIfdiff)); +nameEuler = sprintf('Euler/%.0e', eulerStep); +lwIfdiff = 3; +lwEuler = 3; +colorIfdiff = [0, 0.447, 0.741]; +colorEuler = [0.85, 0.325, 0.098]; +lsIfdiff = '-'; +lsEuler = '--'; + +% Sensitivity +wrt_y = 3; +eulerDisturbH = 1e-6; + + +%% Solve with IFDIFF +datahandle = prepareDatahandleForIntegration(rhs, 'solver', intIfdiff, 'options', intOptions); +% Enable storing of convexification data +data = datahandle.getData(); +data.sliding.storeSlidingInfo = true; +datahandle.setData(data); + +solIfdiff = solveODE(datahandle, tspan, x0, p); + + +%% Solve with Euler +[owndir, ~] = fileparts(mfilename('fullpath')); +if contains(owndir, 'Editor_') + [owndir, ~] = fileparts(matlab.desktop.editor.getActiveFilename); +end + +euler_prefix = sprintf('sol_euler_%.0e', eulerStep); +euler_vname = 'sol_euler'; +euler_fname = fullfile(owndir, [euler_prefix, '.mat']); + +solEuler = loadEulerOrCompute(euler_fname, euler_vname, rhs, tspan, x0, p, eulerStep); + + +%% Plot solution +plot_x_ifdiff = deval(solIfdiff, plot_t); +axSol = plotSol3d([], plot_x_ifdiff, nameIfdiff, lwIfdiff, colorIfdiff, lsIfdiff); + +plot_x_euler = interp1(solEuler.x, solEuler.y', plot_t)'; +plotSol3d(axSol, plot_x_euler, nameEuler, lwEuler, colorEuler, lsEuler); + + +%% Compute Sensitivity with IFDIFF +FDstep = generateFDstep(length(x0), length(p)); +sensFunVDE = generateSensitivityFunction(datahandle, solIfdiff, FDstep, ... + 'method', 'VDE', ... + 'calcGy', true, ... + 'calcGp', false); +sensVDE = sensFunVDE(plot_t); +G = arrayfun(@(s) s.Gy(:, wrt_y), sensVDE, 'UniformOutput', false); +plot_sens_ifdiff = [G{:}]; + + +%% Compute Sensitivity with Euler +euler_prefix_disturb = sprintf('%s_disturb_wrt_y%d', euler_prefix, wrt_y); +euler_fname_disturb = fullfile(owndir, [euler_prefix_disturb, '.mat']); + +x0_disturb = x0; +x0_disturb(wrt_y) = x0_disturb(wrt_y) + eulerDisturbH; +solEulerDisturb = loadEulerOrCompute(euler_fname_disturb, euler_vname, rhs, tspan, x0_disturb, p, eulerStep); +sensEuler = (solEulerDisturb.y - solEuler.y) / eulerDisturbH; +plot_sens_euler = interp1(solEuler.x, sensEuler', plot_t)'; + + +%% Plot Sensitivity +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); +end + + +%% Plot alpha and switching function +data = datahandle.getData(); +t_alpha = data.sliding.convexification.t; +alpha = data.sliding.convexification.alpha; +ax = plotSwitchingFuncOrAlpha([], t_alpha, alpha, 'Alpha', 2, colorIfdiff, 'left', 'Alpha'); + +sw = solIfdiff.switchingFunction{1}; +sw_x = arrayfun(@(t) sw([], t, deval(solIfdiff, t), p), t_alpha); +plotSwitchingFuncOrAlpha(ax, t_alpha, sw_x, 'Switching Function', 2, colorEuler, 'right', 'Switching Function'); + +plotSwitches(ax, solIfdiff.switches) + + +%% Helpers +function sol_euler = loadEulerOrCompute(fname, vname, rhs, tspan, x0, p, step) +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('Saving result to %s for later reuse.\n', fname); +save(fname, vname); +end + +function plotSwitches(ax, switches) +xline(ax, switches, '--', 'LineWidth', 1.5, 'HandleVisibility', 'off'); +end + +function ax = plotSwitchingFuncOrAlpha(ax, t, x, name, lw, color, yyax, ylab) +if isempty(ax) + f = figure; + ax = axes(f); +end +yyaxis(ax, yyax); +hold(ax, 'on'); +plot(ax, t, x, 'DisplayName', name, 'LineWidth', lw, 'Color', color); +hold(ax, 'off'); +grid(ax, 'on'); +xlabel(ax, 't'); +ylabel(ax, ylab); +legend(ax, 'location', 'northeast'); + +yl = ylim(ax); +yl(1) = 0; +ylim(ax, yl); +end + +function ax = plotSens(ax, t, x, name, lw, color, ls) +if isempty(ax) + f = figure; + ax = axes(f); +end +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') +legend(ax, 'location', 'northeast'); +end + +function ax = plotSol3d(ax, x, name, lw, color, ls) +if isempty(ax) + f = figure; + ax = axes(f); +end +hold(ax, 'on'); +plot3(ax, x(3, :), x(2, :), x(1, :), 'DisplayName', name, 'LineWidth', lw, 'Color', color, 'LineStyle', ls); +hold(ax, 'off'); + +view(ax, [97 51]); +box(ax, 'on'); +grid(ax, 'on'); +xlabel(ax, 'Predator'); +ylabel(ax, 'Prey 2'); +zlabel(ax, 'Prey 1'); +legend(ax, 'location', 'northeast'); +end + +function sol = explEuler(rhs, tspan, x0, p, stepsize) +xdim = length(x0); +stepcount = (tspan(end)-tspan(1)) / stepsize; +sfac = 0.001; % store factor +n_out = ceil(stepcount*sfac) + 1; + +Xi = reshape(x0, [], 1); +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); + end + if ~mod(i, stepcount / 100), fprintf('.'); end +end +fprintf('\n'); +T = linspace(tspan(1), tspan(end), n_out); +sol.x = T; +sol.y = X; +end diff --git a/toolbox/examples/predatorpreyFilippov/testpprhs.m b/toolbox/examples/predatorpreyFilippov/testpprhs.m deleted file mode 100644 index 1dc27952..00000000 --- a/toolbox/examples/predatorpreyFilippov/testpprhs.m +++ /dev/null @@ -1,162 +0,0 @@ -%% 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 = 1000; -figure(fignum); clf; hold('on'); -plotit = @plotter; - -% solver selection and configuration -intEuler = @explEuler; -eulerStep = 1e-7; -intIfdiff = @ode45; -intOptions = odeset('reltol', 1e-5, 'abstol', 1e-5, 'MaxStep', 0.5); - -% select what to do -doIfdiff = true; -doEuler = false; % Euler takes very long! false: loads from file, if exists -doErrorplot = true; % compare ifdiff and plain euler - -% name generators -nameIfdiff = @(f) sprintf('ifdiff/%s', func2str(f)); -namePlain = @(f) sprintf('plain %s' , func2str(f)); - - -%% IFDIFF -if doIfdiff - fprintf('Integrating with IFDIFF/%s ...\n', func2str(intIfdiff)) - figure(fignum); - datahandle = prepareDatahandleForIntegration('pprhs', 'solver', func2str(intIfdiff), 'options', intOptions); - th = tic(); - sol_ifdiff = solveODE(datahandle, tspan, x0_1, p); - time_ifdiff = toc(th); fprintf('IFDIFF took %g s\n', time_ifdiff); - X_ifdiff = X_plot; - Y_ifdiff = deval(sol_ifdiff, X_ifdiff); - linewidth = 3.0; - hIFDIFF = plotit(fignum, Y_ifdiff, 'g', nameIfdiff(intIfdiff), linewidth); -end - - -%% EULER Integration -[owndir, ~] = fileparts(mfilename('fullpath')); -euler_fname = fullfile(owndir, sprintf('sol_euler_%.0e.mat', eulerStep)); -if doEuler - 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); - save(euler_fname, "sol_euler"); -elseif isfile(euler_fname) - fprintf('Loading sol_euler from file %s\n', euler_fname); - tmp = load(euler_fname, 'sol_euler'); - sol_euler = tmp.sol_euler; - doEuler = true; -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', namePlain(intEuler), linewidth); -end - -%% ANALYSIS -if doErrorplot - if (doIfdiff && doEuler) - fignum = fignum + 1; errorPlot(fignum, X_plot, Y_ifdiff, Y_euler, nameIfdiff(intIfdiff), namePlain(intEuler)) - else - fprintf('Not enough data to compare. Either IFDIFF or Euler solution missing.\n') - end -end - -% FINITO -return - - -%% HELPERS -function errorPlot(fignum, x1, y1, y2, intname1, intname2) - figure(fignum); clf(fignum); - ydiff = calcDiff(y1, y2); - semilogy(x1, ydiff, 'LineWidth', 1.0); - xlabel('t'); - ylabel('||y||_2'); - title(sprintf('difference %s and %s', intname1, intname2)); - drawnow -end - - -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 - - -function diffnorm = calcDiff(yA, yB) - % Go through x-coordinates of yA, find nearest x-coordinate in yB and compare. - % Algorithm: 1) Go through x of yA in Index i - % 2) Find nearest x in yB in a window centered around yB(:,i) - len = length(yA); - diffnorm = zeros(len, 1); - for i = 1:len - y = yA(:, i); - % find nearest x in yB in a time-window around the timepoint of y - window = 250; - j0 = max(1, i-window); - jf = max(len, j0+window); - jidx = j0:jf; - tmpdiff = yB(:, jidx) - repmat(y, [1, length(jidx)]); - tmpdiff = vecnorm(tmpdiff,2,1); - diffnorm(i) = min(tmpdiff); - diffnorm(i) = diffnorm(i) / norm(y); - end -end - - -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 diff --git a/toolbox/internal/rhsPreprocessing/createDatahandle.m b/toolbox/internal/rhsPreprocessing/createDatahandle.m index 143fe1a1..daa2f331 100644 --- a/toolbox/internal/rhsPreprocessing/createDatahandle.m +++ b/toolbox/internal/rhsPreprocessing/createDatahandle.m @@ -58,5 +58,7 @@ data.caseCtrlif = config.caseCtrlif.default; +data.sliding.storeSlidingInfo = config.storeSlidingInfo; + datahandle = makeClosure(data); end diff --git a/toolbox/internal/solving/initDatahandleFields.m b/toolbox/internal/solving/initDatahandleFields.m index 344e95ae..802aedb8 100644 --- a/toolbox/internal/solving/initDatahandleFields.m +++ b/toolbox/internal/solving/initDatahandleFields.m @@ -16,7 +16,7 @@ function initDatahandleFields(datahandle, tspan, initialvalues, parameters) data.forcedBranching = initDatahandleFields_forcedBranching(); data.getSignature = initDatahandleFields_getSignature(); data.computeSensitivity = initDatahandleFields_computeSensitivity(); -data.sliding = initDatahandleFields_Filippov(); +data.sliding = initDatahandleFields_Filippov(data.sliding); datahandle.setData(data); diff --git a/toolbox/internal/solving/initDatahandleFields_Filippov.m b/toolbox/internal/solving/initDatahandleFields_Filippov.m index 4972f940..7ab85e3f 100644 --- a/toolbox/internal/solving/initDatahandleFields_Filippov.m +++ b/toolbox/internal/solving/initDatahandleFields_Filippov.m @@ -1,18 +1,15 @@ -function sliding = initDatahandleFields_Filippov() +function sliding = initDatahandleFields_Filippov(sliding) %sliding = INITDATAHANDLEFIELDS_FILIPPOV() % %Initialize struct for Filippov mode handling in IFDIFF. % -%OUTPUT: +%INPUT/OUTPUT: % sliding - Struct used for Filippov mode handling in IFDIFF. % struct % %See also INITDATAHANDLEFIELDS. -sliding = extendODE_filippov_regime_cleanup(struct()); - -% if true: convexification parameters and signatures are stored for every integration step -sliding.storeSlidingInfo = false; +sliding = extendODE_filippov_regime_cleanup(sliding); convexification.t = []; convexification.index = []; diff --git a/toolbox/internal/solving/solveODE_setFilippovRHS.m b/toolbox/internal/solving/solveODE_setFilippovRHS.m index 1306e33e..3cf1012a 100644 --- a/toolbox/internal/solving/solveODE_setFilippovRHS.m +++ b/toolbox/internal/solving/solveODE_setFilippovRHS.m @@ -51,8 +51,6 @@ function solveODE_setFilippovRHS(datahandle) % set filippov rhs and store some other info data.sliding.filippov_rhs = filippov_rhs; data.sliding.index = sliding_index; -data.sliding.signature_fplus = signature_fplus; -data.sliding.signature_fminus = signature_fminus; datahandle.setData(data); diff --git a/toolbox/makeConfig.m b/toolbox/makeConfig.m index d235daf1..54deab6f 100644 --- a/toolbox/makeConfig.m +++ b/toolbox/makeConfig.m @@ -30,6 +30,9 @@ % Flag to set debug mode which may display additional information and warnings among other things. config.debugMode = false; +% Flag to enable/disable storing of Filippov convexification parameters for debugging (can be costly!). +config.storeSlidingInfo = false; + % Flag to enable/disable reusing of switching functions that already exist based on the signature of a switch. % Always keep this enabled for best performance. Can be useful to disable for debugging and tests. config.reuseSwitchingFunctions = true;