From 0aca75491e3da86d3930481a73251e7d441443c9 Mon Sep 17 00:00:00 2001 From: pilar Date: Tue, 3 Feb 2026 13:16:34 +0100 Subject: [PATCH 1/2] New instructions (README and webpage) --- README.md | 141 ++++++++++++----------------------------------- docs/index.md | 149 +++++++++++++++++++++++++++++++------------------- 2 files changed, 129 insertions(+), 161 deletions(-) diff --git a/README.md b/README.md index 1c7f27d..2b1bd9a 100644 --- a/README.md +++ b/README.md @@ -1,132 +1,61 @@ -# IFDIFF - A Matlab Toolkit for ODEs with State˗Dependent Switches +# IFDIFF - A MATLAB Toolkit for ODEs with State˗Dependent Switches [![Open in MATLAB Online](https://www.mathworks.com/images/responsive/global/open-in-matlab-online.svg)](https://matlab.mathworks.com/open/github/v1?repo=andreassommer/ifdiff&file=toolbox/doc/GettingStarted.mlx) -The software package IFDIFF deals with the solution and algorithmic generation of sensitivities in -ordinary differential equations with implicit (state-dependent) non-differentiabilites ("switches") -in the right-hand side that is given as Matlab program code with non-differentiable operators such as -`min`, `max`, `abs`, `sign`, as well as `if`-branching. IFDIFF automatically generates only necessary -switching functions, outputs them as Matlab code, and detects switching points accurately up to -machine precision. +The software package IFDIFF comprises: +- **Automatic detection and processing** of state-dependent switching events in ODE IVPs +- **Automatic generation** of only the necessary **switching functions** (exported as MATLAB code) +- Accurate switching point detection up to machine precision +- Detection and handling of **Filippov sliding modes** +- ODEs with state-jumps and model switches +- **Sensitivity generation** for switched ODEs, ODEs with state-jumps and Filippov ODEs -Sensitivities can be generated w.r.t. the initial values and w.r.t. a given parameter set +Sensitivities can be generated w.r.t. both initial values and parameters sets. -See the [IFDIFF page](https://andreassommer.github.io/ifdiff/) for a mathematical introduction with example. +For a mathematical introduction and illustrative examples, see the +[IFDIFF project page](https://andreassommer.github.io/ifdiff/). -The file [Readme_Example.m](./toolbox/examples/Readme_Example.m) contains a self-explaining Matlab script similar to the - contents below. +A compact, self-explanatory MATLAB example is provided in this file. ([`Readme_Example.m`](./toolbox/examples/Readme_Example.m)) -
- - - -# First Run Prerequisites - -Before using IFDIFF, it is mandatory to run the `make_mtreeplus` script once. -After starting Matlab, change to the IFDIFF directory and type `make_mtreeplus`. -This scripts generates a modified copy of Matlab's own parser class `mtree`, on which IFDIFF heavily relies. - -It is advisable to initialize the paths needed for IFDIFF by invoking -`initIFDIFF();` once on every new matlab session.
+# Installation -# Usage - -Make sure that you've followed the [First Run Prerequisites](#first-run-prerequisites). - -The following commands exemplarily show the usage on the canonical example. - -## 1. Preprocessing of Right Hand Side - -The preprocessing analyses and transforms the right hand side function `canonicalExampleRHS.m`. -Preprocessing generates the `datahandle`, the central structure for switching detection and handling. -We set ODE solver and its options as usual. If not set, default values are used. +There are two ways to install IFDIFF. The required [First Run Prerequisites](#first-run-prerequisites) are identical for both methods after installation. - ```matlab - initIFDIFF(); % Initialise the paths for ifdiff -- needed only once per Matlab session - integrator = @ode45; - odeoptions = odeset('AbsTol', 1e-14, 'RelTol', 1e-12); - datahandle = prepareDatahandleForIntegration('canonicalExampleRHS', 'integrator', integrator, 'options', odeoptions); - ``` +## Installation Method A +1. Download the current release file `IFDIFF_Toolbox.mltbx`. +2. Open MATLAB and navigate to the directory containing the file. +3. Right-click the file and select **Open**. +4. For a detailed introduction, open `GettingStarted.mlx`, which is included with the toolbox. -## 2. Integration (Forward Solution) +Note: You can also open the Getting Started guide directly in MATLAB Online and then continue with the [First Run Prerequisites](#first-run-prerequisites). -Define initial values, parameter values, and integration horizon, and call `solveODE` to start the integration. -`solveODE` returns a Matlab `sol` structure, that can be evalated using `deval`. -It is an augmented version of the solution structures returned by Matlab's very own integrator suite -(see https://de.mathworks.com/help/matlab/ref/deval.html#bu7iw_j-sol) - ```matlab - tspan = [0 20]; - initialvalues = [1;0]; - parameters = 5.437; - sol = solveODE(datahandle, tspan, initialvalues, parameters); - T = 0:0.1:20; - X = deval(sol,T); - plot(T,X) -``` +## Installation Method B +1. Navigate to a location of your choice on your PC (e.g. Desktop). +2. Clone the repository `git clone https://github.com/andreassommer/ifdiff.git` +3. Open MATLAB and navigate to the cloned `ifdiff` directory. -## 3. Sensitivity Generation - -The sensitivity generation currently supports 1st order forward sensitivities w.r.t. initial state and parameters. -It is possible to generate sensitivities using external numerical differentiation (flags `END_full`, `END_piecewise`) -or using the variational differential equations (flag `VDE`). - -Usually, method `VDE` delivers vest results in terms of accuracy, as it calculates error-controlled sensitivities. -It uses the variational differential equations on each interval and performes updates at the switching points to -ensure accurate sensitivities at the precalculated forward solution. However, the occuring augmented differential -system might be large size and thus slow to compute. -Required state derivatives of the right hand side are approximated using automated finite differencing. - -The method `END_piecewise` computes the interval sensitivities using external numerical differentiation (finite differencing) -and connects these using the same updates as used in the VDE method. This requires multiple evaluations of interval solutions, -but might be faster on larger systems. - -When using `END_full`, external numerical differention is used on multiple full horizon solutions. The individual -trajectories are calculated with switching point detection. Thus, no updates between switches are required. -In general, this approach is less accurate and slower than `END_piecewise`, but might be the a good choice -for highly instable ODEs. - -1. Choose step sizes for finite differencing (also used in method `VDE` for generating state derivatives of the rhs). - - ```matlab - dim_y = size(sol.y,1); // state dimension - dim_p = length(parameters); // number of parameters - FDstep = generateFDstep(dim_y,dim_p); - ``` - - The `generateDFstep` function accepts several options influencing e.g. step length. - See the documentation or the file for more information +
+# First Run Prerequisites -2. Build the sensitivity function. In this example, the `END_piecewise` method is chosen. +Before using IFDIFF, the script `make_mtreeplus` must be executed once. (In `ifdiff/toolbox` for installation Method B.) +This script generates a modified copy of MATLAB's internal parser class `mtree`, which is required by IFDIFF. - ```matlab - sensitivity_function = generateSensitivityFunction(datahandle, sol, FDstep, 'method', 'END_piecewise'); - ``` +In each new MATLAB session, initialize the required paths by calling: - The following table lists several name-value pairs that can be used to configure `generateSensitivityFunction` +`initIFDIFF();` - | Parameters | Possible values | Defaults | - | ----------------------- | ----------------------------------------------------------------------------------- | ------------------------------------- | - | calcGy | true/false - flag indicating to calculate state sensitivities | true | - | calcGp | true/false - flag indicating to calculate parameter sensitivities | true | - | Gmatrices_intermediate | true/false - flag indicating to store update matrics | false | - | save_intermediates | true/false - flag indicating to store intermediate calculations | true | - | integrator | Function handle for ODE solver in Matlab (e.g. ode45) | Integrator used by ifdiff | - | integrator_options | Options struct generated for ODE solver | Integrator options used by ifdiff | - | method | String with VDE/END_piecewise/END_full | VDE | - | directions_y | Matrix containing directions for directional derivatives w.r.t initial values. | Identity matrix with dimension n_y | - | directions_p | Matrix containing directions for directional derivatives w.r.t parameters. | Identity matrix with dimension n_p | +This command can be executed either from the toolbox directory or from the cloned repository. +
-3. Evaluate the sensitivity function at specific times. +# Usage - ```matlab - t = 0:0.1:20; - sensitivities = sensitivity_function(t); - ``` +Ensure that the [First Run Prerequisites](#first-run-prerequisites) have been completed. +For a step-by-step usage example and background explanation, see the [IFDIFF project page](https://andreassommer.github.io/ifdiff/). diff --git a/docs/index.md b/docs/index.md index 29cccf5..37b1366 100644 --- a/docs/index.md +++ b/docs/index.md @@ -1,21 +1,17 @@ -# IFDIFF - A Matlab Toolkit for ODEs with State˗Dependent Switches +# IFDIFF - A MATLAB Toolkit for ODEs with State˗Dependent Switches -The software package IFDIFF deals with the solution and algorithmic generation of sensitivities -in ordinary differential equations with implicit (state-dependent) non-differentiabilites ("switches") -in the right-hand side that is given as Matlab program code with non-differentiable operators such as -`min`, `max`, `abs`, `sign`, as well as `if`-branching. IFDIFF automatically generates only necessary -switching functions, outputs them as Matlab code, and detects switching points accurately up to machine -precision. -IFDIFF handles multidimensional state and parameter vectors and can be transparently used in existing code, -as it generates `sol` solution structures compatible to the Matlab ode solvers. No modifications to the -right hand side or manual processing are necessary. A single preparation call is sufficient. +The software package **IFDIFF** provides algorithms for the accurate numerical solution and sensitivity analysis of ordinary differential equations (ODEs) with implicit, state-dependent non-differentiabilities ("switches") in the right-hand side as well as Filippov ODEs and ODEs with state-jumps. -In the PErFDiff project, the existing implementation is to be further developed and extended for Filippov systems. +The right-hand side is supplied as standard MATLAB program code and may contain non-differentiable operators such as `min`, `max`, `abs`, `sign` or `if`-branching. +IFDIFF automatically detects and processes these switches, generates only the necessary switching functions (exported as MATLAB code), and determines switching points accurately up to machine precision. -
-
+IFDIFF handles multidimensional state and parameter vectors and produces solution structures compatible with MATLAB's built-in ODE solvers, allowing transparent use of functions such as `deval` without modification of existing code. +A single preparation call is sufficient to enable switching-aware integration. +Within the **PErFDiff** project, the existing IFDIFF implementation is being further developed and extended to cover Filippov systems. + +
# Canonical Example Consider the following "canonical example" for a switched ODE system: @@ -27,7 +23,7 @@ $$ with $$ - f_1(t,x,p) = 0.01 * t^2 + x_2^3 + f_1(t,x,p) = 0.01 \cdot t^2 + x_2^3 \qquad f_2(t,x,p) = \begin{cases} 0 ~~~if~~ x_1 < p \\ 5 ~~~if~~ p \leq x_1 < p+0.5 \\ 0 ~~~if~~ x_1 \geq p+0.5 \end{cases} $$ @@ -58,12 +54,11 @@ and a second "kink" once the first component reaches `p+0.5`. Let's see what happens if we do not consider appropriate switching handling. -
-
+
# Naive Integration with ode45 fails unnoticed -We initialize the variables and start the integration using Matlab's default integrator `ode45` (explicit Runge-Kutta 4(5)-solver), +We initialize the variables and start the integration using MATLAB's default integrator `ode45` (explicit Runge-Kutta 4(5)-solver), without caring for the non-differentiable `if` statements. ```matlab tspan = [0 20]; % time horizon @@ -83,8 +78,7 @@ Contrary to wide-spread beliefs, tightening the integration tolerances is not a Without any warning or error, the __naive approach with tight tolerances leads to the same (wrong) result__. Again, no kinks visible. -
-
+
# Reliable integration with IFDIFF @@ -107,36 +101,89 @@ are exact up to integration tolerance. Using, e.g., absolute and relative tolera first switching times is computed an error less than $10^{-14}$, and the second with an error less than $10^{-11}$. Note that the `sol` structure returned by `solveODE` is an augmented version of the solution structures returned -by Matlab's very own integrators (see [Matlab documentation](https://de.mathworks.com/help/matlab/ref/deval.html#bu7iw_j-sol)), +by MATLAB's very own integrators (see [MATLAB documentation](https://de.mathworks.com/help/matlab/ref/deval.html#bu7iw_j-sol)), and can thus be evaluated using `deval`. That means, IFDIFF can be used transparently within existing code! -
+
## Remarks -While in many (not all) cases, forcing a small integration step size (e.g. using `odeset('MaxStep',1e-2)` in the above example), -may lead to approximative solutions, this is not a remedy at all. +While in some cases forcing a small integration step size (for example by using `odeset('MaxStep',1e-2)` in the example above) may produce a visually plausible solution, this approach is not a reliable remedy. + +Firstly, restricting the maximum step size contradicts the fundamental idea of adaptive step-size integrators. It globally limits the step size, even in regions where no switching occurs. +In the example above, choosing a parameter value of `p = 100` would eliminate all switches within the integration horizon, yet the step size would still be unnecessarily constrained. + +Secondly, when non-differentiabilities are ignored, there is no principled way to choose an appropriate maximum step size. The integrator provides neither warnings nor error estimates indicating that switching events have been missed. + +Thirdly, and as demonstrated in the following section, this approach completely fails to produce meaningful sensitivities. + +
+ +# Sensitivity Generation + +Correct state and parameter sensitivities are crucial for derivative-based optimization and control. +While restricting the maximum integrator step size (e.g. to `0.1`) may allow the naive approach to reproduce the forward trajectory, the resulting sensitivities are completely unreliable and generally unusable. + +IFDIFF currently supports first-order forward sensitivities w.r.t. both the initial state and model parameters. Sensitivities can be computed using either external numerical differentiation or variational differential equations: + +- **Variational Differential Equations** (`VDE`) +This method typically provides the highest accuracy, as it computes error-controlled sensitivities alongside the forward solution. The variational equations are integrated on each interval between switching events, with appropriate updates applied at the detected switching points to ensure consistency with the forward trajectory. +The required derivatives of the right-hand side are approximated using automated finite differencing. +Due to the augmented system size, this method may become computationally expensive for large-scale problems. + +- **External Numerical Differentiation, Piecewise** (`END_piecewise`) +In this approach, sensitivities are computed via finite differencing on each interval between switching events and then connected using the same update rules as in the VDE method. +Although this requires multiple forward integrations per interval, it can be more efficient than VDE for larger systems while still providing accurate sensitivities. -First, restricting the maximum step size is completely against the idea of adaptive step size integrators, as it reduces the maximum -step size everywhere - even if there is no switch at all (consider, in the above example, a parameter value of `p=100`; then there is -no switch in the specified integration horizon). +- **External Numerical Differentiation, Full Horizon** (`END_full`) +Here, finite differencing is applied to multiple full-horizon solutions, each computed with switching point detection. Since sensitivities are not propagated across switching points, no intermediate updates are required. +This method is generally less accurate and we do not recommend it, but may be preferable for highly unstable ODEs where interval-based propagation becomes unreliable. -Second, using the naive approach (i.e. ignoring non-differentiabilities completely), we don't know how to choose ths maximum step size -appropriately, as we don't get any error messages or warnings at all. +The following table lists several name-value pairs that can be used to configure `generateSensitivityFunction` -Third, and this we will show in the next section, sensitivites cannot be computed at all by this approach. + | Parameters | Possible values | Defaults | + | ----------------------- | ----------------------------------------------------------------------------------- | ------------------------------------- | + | calcGy | true/false - flag indicating to calculate state sensitivities | true | + | calcGp | true/false - flag indicating to calculate parameter sensitivities | true | + | Gmatrices_intermediate | true/false - flag indicating to store update matrics | false | + | save_intermediates | true/false - flag indicating to store intermediate calculations | true | + | integrator | Function handle for ODE solver in MATLAB (e.g. ode45) | Integrator used by ifdiff | + | integrator_options | Options struct generated for ODE solver | Integrator options used by ifdiff | + | method | String with VDE/END_piecewise/END_full | VDE | + | directions_y | Matrix containing directions for directional derivatives w.r.t initial values. | Identity matrix with dimension n_y | + | directions_p | Matrix containing directions for directional derivatives w.r.t parameters. | Identity matrix with dimension n_p | -
-
+
-# Computing Sensitivities +## Example sensitivity generation for the Canonical Example -Having correct state and parameters sensitivities, i.e. the derivative of the trajectory w.r.t. changes in the initial values or parameters, -is cruicial for derivative-based optimization. -By restricting the maximum integrator step size to $0.1$, also the naive approach is capable of reproducing the forward trajectory. -However, the sensitivities generated by the naive approach are completely useless. +1. Choose step sizes for finite differencing (also used in method `VDE` for generating state derivatives of the rhs). -The following picture shows the trajectories of the state sensitivities, $G_{y,ij}(t,t_0) := \frac{d y_i}{d x_{0,j}}(t)$, i.e. + ```matlab + dim_y = size(sol.y,1); // state dimension + dim_p = length(parameters); // number of parameters + FDstep = generateFDstep(dim_y,dim_p); + ``` + + The `generateDFstep` function accepts several options influencing e.g. step length. + See the documentation or the file for more information + + +2. Build the sensitivity function. In this example, the `END_piecewise` method is chosen. + + ```matlab + sensitivity_function = generateSensitivityFunction(datahandle, sol, FDstep, 'method', 'END_piecewise'); + +3. Evaluate the sensitivity function at specific times. + + ```matlab + t = 0:0.1:20; + sensitivities = sensitivity_function(t); + ``` + +
+ +The following picture shows the trajectories of the state sensitivities for the example aboves, $G_{y,ij}(t,t_0) := \frac{d y_i}{d x_{0,j}}(t)$, i.e. $G_{y,ij}$ denotes the sensitivity of the $i$-th solution component w.r.t. the $j$-th component of the initial value. | Comparison: Sensitivities w.r.t. initial state | | @@ -144,30 +191,22 @@ $G_{y,ij}$ denotes the sensitivity of the $i$-th solution component w.r.t. the $ | ![Sensitivities with naive approach](./canonex_sensitivity_naive.png) | ![Sensitivities with ifdiff](./canonex_sensitivity_ifdiff.png) | | __Naive approach__: Sensitivities are useless | __IFDIFF__: Correct and accurate sensitivities | -
+
While IFDIFF produces correct and accurate sensitivities, the sensitivities using the naive approach, again without any warning or hint, are not only inaccurate (notice the scales!), but plain wrong. Details of how to compute sensitivities with IFDIFF are given in the project's [README.md](https://github.com/andreassommer/ifdiff/blob/public/README.md). -
-
+
# IFDIFF is simple to use! -Correct treatment of switched systems requires elaborate formulation of switching functions and tailored integrators, -placing high mathematical demands on modelers. Even small model changes often imply considerable reformulation effort. -Furthermore: $n$ switches generate up to $2^n$ possible program flows and switching functions, rendering a-priori -formulations not feasible already in medium-sized models. - -IFDIFF programmatically handles switching events, auto-generating only required switching functions. -It determines switching times up to machine precision, and ensures accurate simulation and sensitivity results. -Transparently extending the Matlab integrators (ode45, ode15s, etc.), IFDIFF is applicable to existing code with -state- and parameter-dependent conditionals, thus enabling fast prototyping and relieving modelers of -mathematical-technical effort. - -Calculation times of the naive approach are a little lower than the ones of IFDIFF. -But it can generate arbitrarily wrong results without any notice (see the example above!). -IFDIFF provides correct integration results, correct first-order sensitivities and information -about the switching structure of your model. +Accurate simulation of switched systems typically requires explicit formulation of switching functions and specialized integrators, placing high mathematical demands on the modeler. +Even minor model changes can require substantial reformulation, and with $n$ switches up to $2^n$ possible execution paths may arise, making a priori formulations infeasible for medium-sized systems. + +IFDIFF handles switching events programmatically, automatically generating only the required switching functions and detecting switching times up to machine precision. +By extending MATLAB's standard ODE integrators (`ode45`, `ode15s`, etc.), IFDIFF can be applied directly to existing code containing state- or parameter-dependent conditionals, enabling fast prototyping without manual reformulation. + +While naive integration may appear slightly faster, it can silently produce incorrect trajectories and meaningless sensitivities. +IFDIFF instead provides reliable solutions, accurate first-order sensitivities, and explicit information about the switching structure of the model. \ No newline at end of file From 95f4e6543270bf697de54a3c8f707c59f965a01a Mon Sep 17 00:00:00 2001 From: pilar Date: Tue, 10 Feb 2026 11:34:32 +0100 Subject: [PATCH 2/2] Reworked instructions --- README.md | 29 +++++------- docs/index.md | 125 +++++++++++++++++++++++++------------------------- 2 files changed, 73 insertions(+), 81 deletions(-) diff --git a/README.md b/README.md index 2b1bd9a..611212b 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ The software package IFDIFF comprises: - **Automatic detection and processing** of state-dependent switching events in ODE IVPs - **Automatic generation** of only the necessary **switching functions** (exported as MATLAB code) - Accurate switching point detection up to machine precision -- Detection and handling of **Filippov sliding modes** +- Detection and handling of **Filippov sliding mode** - ODEs with state-jumps and model switches - **Sensitivity generation** for switched ODEs, ODEs with state-jumps and Filippov ODEs @@ -22,19 +22,19 @@ A compact, self-explanatory MATLAB example is provided in this file. ([`Readme_E # Installation -There are two ways to install IFDIFF. The required [First Run Prerequisites](#first-run-prerequisites) are identical for both methods after installation. +There are two ways to install IFDIFF. -## Installation Method A +## Installation Method (Recommended) 1. Download the current release file `IFDIFF_Toolbox.mltbx`. 2. Open MATLAB and navigate to the directory containing the file. 3. Right-click the file and select **Open**. 4. For a detailed introduction, open `GettingStarted.mlx`, which is included with the toolbox. -Note: You can also open the Getting Started guide directly in MATLAB Online and then continue with the [First Run Prerequisites](#first-run-prerequisites). +Note: You can also open the Getting Started guide directly in +[MATLAB Online](https://matlab.mathworks.com/open/github/v1?repo=andreassommer/ifdiff&file=toolbox/doc/GettingStarted.mlx) and then continue with the [First Run Prerequisites](#first-run-prerequisites). - -## Installation Method B +## Installation Method (Alternative) 1. Navigate to a location of your choice on your PC (e.g. Desktop). 2. Clone the repository `git clone https://github.com/andreassommer/ifdiff.git` @@ -42,20 +42,13 @@ Note: You can also open the Getting Started guide directly in MATLAB Online and
-# First Run Prerequisites - -Before using IFDIFF, the script `make_mtreeplus` must be executed once. (In `ifdiff/toolbox` for installation Method B.) -This script generates a modified copy of MATLAB's internal parser class `mtree`, which is required by IFDIFF. - -In each new MATLAB session, initialize the required paths by calling: - -`initIFDIFF();` +# Usage -This command can be executed either from the toolbox directory or from the cloned repository. +Onnce every new MATLAB session run -
+ `initIFDIFF();` -# Usage +to initialize the required paths. Then IFDIFF is ready to use. -Ensure that the [First Run Prerequisites](#first-run-prerequisites) have been completed. +The first time the command is executed, a copy of MATLAB's internal parser class `mtree` is created on which IFDIFF heavily relies. For a step-by-step usage example and background explanation, see the [IFDIFF project page](https://andreassommer.github.io/ifdiff/). diff --git a/docs/index.md b/docs/index.md index 37b1366..696ace5 100644 --- a/docs/index.md +++ b/docs/index.md @@ -1,7 +1,7 @@ # IFDIFF - A MATLAB Toolkit for ODEs with State˗Dependent Switches -The software package **IFDIFF** provides algorithms for the accurate numerical solution and sensitivity analysis of ordinary differential equations (ODEs) with implicit, state-dependent non-differentiabilities ("switches") in the right-hand side as well as Filippov ODEs and ODEs with state-jumps. +The software package **IFDIFF** provides tools for the accurate numerical solution and sensitivity analysis of ordinary differential equations (ODEs) with state-dependent (implicit) non-differentiabilities ("switches") in the right-hand side as well as Filippov Sliding mode and ODEs with state-jumps. The right-hand side is supplied as standard MATLAB program code and may contain non-differentiable operators such as `min`, `max`, `abs`, `sign` or `if`-branching. IFDIFF automatically detects and processes these switches, generates only the necessary switching functions (exported as MATLAB code), and determines switching points accurately up to machine precision. @@ -9,26 +9,25 @@ IFDIFF automatically detects and processes these switches, generates only the ne IFDIFF handles multidimensional state and parameter vectors and produces solution structures compatible with MATLAB's built-in ODE solvers, allowing transparent use of functions such as `deval` without modification of existing code. A single preparation call is sufficient to enable switching-aware integration. -Within the **PErFDiff** project, the existing IFDIFF implementation is being further developed and extended to cover Filippov systems.
# Canonical Example Consider the following "canonical example" for a switched ODE system: -$$ +$${ \dot x = f(t,x,p) = \binom{f_1(t,x,p)}{f_2(t,x,p)} -$$ +}$$ with -$$ +$${ f_1(t,x,p) = 0.01 \cdot t^2 + x_2^3 \qquad f_2(t,x,p) = \begin{cases} 0 ~~~if~~ x_1 < p \\ 5 ~~~if~~ p \leq x_1 < p+0.5 \\ 0 ~~~if~~ x_1 \geq p+0.5 \end{cases} -$$ +}$$ -and with initial value $ x(0) = (1,0)^T $, parameter $ p = 5.437 $, over time span $ t \in [0,20] $. +and with initial value ${ x(0) = (1,0)^T }$, parameter ${ p = 5.437 }$, over time span ${ t \in [0,20] }$. This switched ODE system translates straightforward into the following matlab program: @@ -60,19 +59,20 @@ Let's see what happens if we do not consider appropriate switching handling. We initialize the variables and start the integration using MATLAB's default integrator `ode45` (explicit Runge-Kutta 4(5)-solver), without caring for the non-differentiable `if` statements. + ```matlab - tspan = [0 20]; % time horizon - x0 = [1;0]; % initial values - p = 5.437; % parameter values - solX = ode45(@(t,x) canonicalExampleRHS(t,x,p), tspan, x0) - T = 0:0.1:20; X = deval(solX2,T); plot(T,X); legend('x_1','x_2'); +tspan = [0 20]; % time horizon +x0 = [1;0]; % initial values +p = 5.437; % parameter values +solX = ode45(@(t,x) canonicalExampleRHS(t,x,p), tspan, x0) +T = 0:0.1:20; X = deval(solX2,T); plot(T,X); legend('x_1','x_2'); ``` ![Canonical Example with naive ode45](./canonex_naive.png) Contrary to wide-spread beliefs, tightening the integration tolerances is not a remedy! ```matlab - odeopts = odeset('AbsTol', 1e-20, 'RelTol', 1e-14); - solX2 = ode45(@(t,x) canonicalExampleRHS(t,x,p), tspan, x0) +odeopts = odeset('AbsTol', 1e-20, 'RelTol', 1e-14); +solX2 = ode45(@(t,x) canonicalExampleRHS(t,x,p), tspan, x0) ``` ![Canonical Example with naive ode45](./canonex_naive_highaccuracy.png) @@ -85,13 +85,16 @@ Without any warning or error, the __naive approach with tight tolerances leads t With the switching point detection in IFDIFF, after a single call to a preparation routine, integration is just as simple as before: ```matlab - initIFDIFF(); % initialise IFDIFF (only once) - tspan = [0 20]; x0 = [1;0]; p = 5.437; % set time horizon, initial value, parameter - integrator = @ode45; % choose integrator - odeoptions = odeset('AbsTol', 1e-5, 'RelTol', 1e-3); % set integrator options, here: low accuracy - datahandle = prepareDatahandleForIntegration('canonicalExampleRHS', 'integrator', func2str(integrator), 'options', odeoptions); - sol = solveODE(datahandle, tspan, x0, p); - T = 0:0.1:20; X = deval(solX2,T); plot(T,X); legend('x_1','x_2'); +initIFDIFF(); % initialise IFDIFF (only once) +tspan = [0 20]; x0 = [1;0]; p = 5.437; % set time horizon, initial value, parameter +integrator = @ode45; % choose integrator +odeoptions = odeset('AbsTol', 1e-5, 'RelTol', 1e-3); % set integrator options, here: low accuracy + +datahandle = prepareDatahandleForIntegration('canonicalExampleRHS', ... + 'integrator', integrator, ... + 'options', odeoptions); +sol = solveODE(datahandle, tspan, x0, p); +T = 0:0.1:20; X = deval(solX2,T); plot(T,X); legend('x_1','x_2'); ``` and __IFDIFF delivers the correct result__: ![Canonical Example with ifdiff](./canonex_ifdiff.png) @@ -135,55 +138,39 @@ Due to the augmented system size, this method may become computationally expensi In this approach, sensitivities are computed via finite differencing on each interval between switching events and then connected using the same update rules as in the VDE method. Although this requires multiple forward integrations per interval, it can be more efficient than VDE for larger systems while still providing accurate sensitivities. -- **External Numerical Differentiation, Full Horizon** (`END_full`) +- **External Numerical Differentiation, Full Horizon** (`END_full`) **NOT Recommended** Here, finite differencing is applied to multiple full-horizon solutions, each computed with switching point detection. Since sensitivities are not propagated across switching points, no intermediate updates are required. -This method is generally less accurate and we do not recommend it, but may be preferable for highly unstable ODEs where interval-based propagation becomes unreliable. - -The following table lists several name-value pairs that can be used to configure `generateSensitivityFunction` - - | Parameters | Possible values | Defaults | - | ----------------------- | ----------------------------------------------------------------------------------- | ------------------------------------- | - | calcGy | true/false - flag indicating to calculate state sensitivities | true | - | calcGp | true/false - flag indicating to calculate parameter sensitivities | true | - | Gmatrices_intermediate | true/false - flag indicating to store update matrics | false | - | save_intermediates | true/false - flag indicating to store intermediate calculations | true | - | integrator | Function handle for ODE solver in MATLAB (e.g. ode45) | Integrator used by ifdiff | - | integrator_options | Options struct generated for ODE solver | Integrator options used by ifdiff | - | method | String with VDE/END_piecewise/END_full | VDE | - | directions_y | Matrix containing directions for directional derivatives w.r.t initial values. | Identity matrix with dimension n_y | - | directions_p | Matrix containing directions for directional derivatives w.r.t parameters. | Identity matrix with dimension n_p | - +This method is generally not accurate and we do not recommend it, but may be slower for a certain class of problems. +
## Example sensitivity generation for the Canonical Example 1. Choose step sizes for finite differencing (also used in method `VDE` for generating state derivatives of the rhs). - ```matlab - dim_y = size(sol.y,1); // state dimension - dim_p = length(parameters); // number of parameters - FDstep = generateFDstep(dim_y,dim_p); - ``` + ```matlab + dim_y = size(sol.y,1); + dim_p = length(parameters); + FDstep = generateFDstep(dim_y,dim_p); + ``` The `generateDFstep` function accepts several options influencing e.g. step length. - See the documentation or the file for more information + See the documentation below for more information. 2. Build the sensitivity function. In this example, the `END_piecewise` method is chosen. - ```matlab - sensitivity_function = generateSensitivityFunction(datahandle, sol, FDstep, 'method', 'END_piecewise'); - + ```matlab + sensitivity_function = generateSensitivityFunction(datahandle, sol, FDstep, 'method', 'END_piecewise'); + ``` 3. Evaluate the sensitivity function at specific times. - ```matlab - t = 0:0.1:20; - sensitivities = sensitivity_function(t); - ``` - -
+ ```matlab + t = 0:0.1:20; + sensitivities = sensitivity_function(t); + `` -The following picture shows the trajectories of the state sensitivities for the example aboves, $G_{y,ij}(t,t_0) := \frac{d y_i}{d x_{0,j}}(t)$, i.e. +The following picture shows the trajectories of the state sensitivities for the example aboves, $G_{y,ij}(t,t_0) := \frac{d y_i}{d y_{0,j}}(t)$, i.e. $G_{y,ij}$ denotes the sensitivity of the $i$-th solution component w.r.t. the $j$-th component of the initial value. | Comparison: Sensitivities w.r.t. initial state | | @@ -191,22 +178,34 @@ $G_{y,ij}$ denotes the sensitivity of the $i$-th solution component w.r.t. the $ | ![Sensitivities with naive approach](./canonex_sensitivity_naive.png) | ![Sensitivities with ifdiff](./canonex_sensitivity_ifdiff.png) | | __Naive approach__: Sensitivities are useless | __IFDIFF__: Correct and accurate sensitivities | -
- While IFDIFF produces correct and accurate sensitivities, the sensitivities using the naive approach, again without any warning or hint, are not only inaccurate (notice the scales!), but plain wrong. -Details of how to compute sensitivities with IFDIFF are given in the project's [README.md](https://github.com/andreassommer/ifdiff/blob/public/README.md). +
+ +The following table lists several name-value pairs that can be used to configure `generateSensitivityFunction` + + | Parameters | Possible values | Defaults | + | ----------------------- | ----------------------------------------------------------------------------------- | ------------------------------------- | + | calcGy | true/false - flag indicating to calculate state sensitivities | true | + | calcGp | true/false - flag indicating to calculate parameter sensitivities | true | + | Gmatrices_intermediate | true/false - flag indicating to store update matrics | false | + | save_intermediates | true/false - flag indicating to store intermediate calculations | true | + | integrator | Function handle for ODE solver in MATLAB (e.g. ode45) | Integrator used by ifdiff | + | integrator_options | Options struct generated for ODE solver | Integrator options used by ifdiff | + | method | String with VDE/END_piecewise/END_full | VDE | + | directions_y | Matrix containing directions for directional derivatives w.r.t initial values. | Identity matrix with dimension n_y | + | directions_p | Matrix containing directions for directional derivatives w.r.t parameters. | Identity matrix with dimension n_p |
+ # IFDIFF is simple to use! -Accurate simulation of switched systems typically requires explicit formulation of switching functions and specialized integrators, placing high mathematical demands on the modeler. -Even minor model changes can require substantial reformulation, and with $n$ switches up to $2^n$ possible execution paths may arise, making a priori formulations infeasible for medium-sized systems. +Correct treatment of switched systems requires elaborate formulation of switching functions and tailored integrators, placing high mathematical demands on modelers. +Even small model changes often imply considerable reformulation effort. Furthermore: ${n}$ switches generate up to ${ 2^n}$ possible program flows and switching functions, rendering a-priori formulations not feasible already in medium-sized models. -IFDIFF handles switching events programmatically, automatically generating only the required switching functions and detecting switching times up to machine precision. -By extending MATLAB's standard ODE integrators (`ode45`, `ode15s`, etc.), IFDIFF can be applied directly to existing code containing state- or parameter-dependent conditionals, enabling fast prototyping without manual reformulation. +IFDIFF programmatically handles switching events, auto-generating only required switching functions. It determines switching times up to machine precision, and ensures accurate simulation and sensitivity results. Transparently extending the Matlab integrators (ode45, ode15s, etc.), IFDIFF is applicable to existing code with state- and parameter-dependent conditionals, thus enabling fast prototyping and relieving modelers of mathematical-technical effort. -While naive integration may appear slightly faster, it can silently produce incorrect trajectories and meaningless sensitivities. -IFDIFF instead provides reliable solutions, accurate first-order sensitivities, and explicit information about the switching structure of the model. \ No newline at end of file +Calculation times of the naive approach are a little lower than the ones of IFDIFF. +But it can generate arbitrarily wrong results without any notice (see the example above!). IFDIFF provides correct integration results, correct first-order sensitivities and information about the switching structure of your model. \ No newline at end of file