Search

We consider piecewise-smooth dynamical systems, i.e., systems of ordinary differential equations switching between different sets of equations on distinct domains, separated by hyper-surfaces. As is well-known, when the solution approaches a discontinuity manifold, a classical solution may cease to exist. For this reason, starting with the pioneering work of Filippov, a concept of weak solution (also known as sliding mode) has been introduced and studied. Nowadays, the solution of piecewise-smooth dynamical systems in and close to discontinuity manifolds is well understood, if the manifold consists locally of a single discontinuity hyper-surface or of the intersection of two discontinuity hyper-surfaces. The present work presents partial results on the solution in and close to discontinuity manifolds of codimension 3 and higher.

In the wake of the 2020 COVID-19 epidemic, much work has been performed on the development of mathematical models for the simulation of the epidemic, and of disease models generally. Most works follow the susceptible-infected-removed (SIR) compartmental framework, modeling the epidemic with a system of ordinary differential equations. Alternative formulations using a partial differential equation (PDE) to incorporate both spatial and temporal resolution have also been introduced, with their numerical results showing potentially powerful descriptive and predictive capacity. In the present work, we introduce a new variation to such models by using delay differential equations (DDEs). The dynamics of many infectious diseases, including COVID-19, exhibit delays due to incubation periods and related phenomena. Accordingly, DDE models allow for a natural representation of the problem dynamics, in addition to offering advantages in terms of computational time and modeling, as they eliminate the need for additional, difficult-to-estimate, compartments (such as exposed individuals) to incorporate time delays. In the present work, we introduce a DDE epidemic model in both an ordinary- and partial differential equation framework. We present a series of mathematical results assessing the stability of the formulation. We then perform several numerical experiments, validating both the mathematical results and establishing model’s ability to reproduce measured data on realistic problems.

We propose and study an algorithm for computing a nearest passive system to a given nonpassive linear time-invariant system (with much freedom in the choice of the metric defining “nearest,” which ...

Objective COVID-19 data released by public health authorities is subject to inherent time delays. Such delays have many causes, including delays in data reporting and the natural incubation period of the disease. We develop and introduce a numerical procedure to recover the distribution of these delays from data. Methods We extend a previously-introduced compartmental model with a nonlinear, distributed-delay term with a general distribution, obtaining an integrodifferential equation. We show this model can be approximated by a weighted-sum of constant time-delay terms, yielding a linear problem for the distribution weights. Standard optimization can then be used to recover the weights, approximating the distribution of the time delays. We demonstrate the viability of the approach against data from Italy and Austria. Results We find that the delay-distributions for both Italy and Austria follow a Gaussian-like profile, with a mean of around 11 to 14 days. However, we note that the delay does not appear constant across all data types, with infection, recovery, and mortality data showing slightly different trends, suggesting the presence of independent delays in each of these processes. We also found that the recovered delay-distribution is not sensitive to the discretization resolution. Conclusions These results establish the validity of the introduced procedure for the identification of time-delays in COVID-19 data. Our methods are not limited to COVID-19, and may be applied to other types of epidemiological data, or indeed any dynamical system with time-delay effects.

In this article, a new method is proposed to approximate the rightmost eigenpair of certain matrix-valued linear operators, in a low-rank setting. First, we introduce a suitable ordinary differential equation, whose solution allows us to approximate the rightmost eigenpair of the linear operator. After analyzing the behaviour of its solution on the whole space, we project the ODE on a low-rank manifold of prescribed rank and correspondingly analyze the behaviour of its solutions. For a general linear operator we prove that—under generic assumptions—the solution of the ODE converges globally to its leading eigenmatrix. The analysis of the projected operator is more subtle due to its nonlinearity; when ca is self-adjoint, we are able to prove that the associated low-rank ODE converges (at least locally) to its rightmost eigenmatrix in the low-rank manifold, a property which appears to hold also in the more general case. Two explicit numerical methods are proposed, the second being an adaptation of the projector splitting integrator proposed recently by Lubich and Oseledets. The numerical experiments show that the method is effective and competitive.

We consider the problem of detecting the common dynamic among several observed signals. It has been shown in (Markovsky et al., 2019) that the problem is equivalent to a generalization of the classical Hankel low-rank approximation to the case of multiple rank constraints. We propose an optimization method based on the integration of ordinary differential equations describing a descent dynamic for a suitable functional to be minimized. We show how the proposed algorithm improves the numerical solutions computed by existing subspace methods which solve the same problem.

In this article, we consider the problems (unsolved in the literature) of computing the nearest normal matrix X to a given non-normal matrix A, under certain constraints, that are (i) if A is real, we impose that also X is real; (ii) if A has known entries on a given sparsity pattern Ω and unknown/uncertain entries otherwise, we impose to X the constraint xij = aij for all entries (i,j) in the pattern Ω. As far as we know, there do not exist in the literature specific algorithms aiming to solve these problems. For the case in which all entries of A can be modified, there exists an algorithm by Ruhe, which is able to compute the closest normal matrix. However, if A is real, the closest computed matrix by Ruhe’s algorithm might be complex, which motivates the development of a different algorithm preserving reality. Normality is characterized in a very large number of ways; in this article, we consider the property that the square of the Frobenius norm of a normal matrix is equal to the sum of the squares of the moduli of its eigenvalues. This characterization allows us to formulate as equivalent problem the minimization of a functional of an unknown matrix, which should be normal, fulfill the required constraints, and have minimal distance from the given matrix A.

A novel algorithmic approach is presented for the problem of partitioning a connected undirected weighted graph under constraints such as cardinality or membership requirements or must-link and can...

This article considers the numerical treatment of piecewise-smooth dynamical systems. Classical solutions as well as sliding modes up to codimension-2 are treated. An algorithm is presented that, in the case of non-uniqueness, selects a solution that is the formal limit solution of a regularized problem. The numerical solution of a regularized differential equation, which creates stiffness and often also high oscillations, is avoided.

We generalize ideas in the recent literature and develop new ones in order to propose a general class of contour integral methods for linear convection–diffusion PDEs and in particular for those arising in finance. These methods aim to provide a numerical approximation of the solution by computing its inverse Laplace transform. The choice of the integration contour is determined by the computation of a few suitably weighted pseudo-spectral level sets of the leading operator of the equation. Parabolic and hyperbolic profiles proposed in the literature are investigated and compared to the elliptic contour originally proposed by Guglielmi, López-Fernández and Nino 2020, see Guglielmi et al. (Math Comput 89:1161–1191, 2020). In summary, the article provides a comparison among three different integration profiles; proposes a new fast pseudospectral roaming method; optimizes the selection of time windows on which one may arbitrarily approximate the solution by no extra computational cost with respect to the case of a fixed time instant; focuses extensively on computational aspects and it is the reference of the MATLAB code [20], where all algorithms described here are implemented.

In this article, we propose a new method for low-rank completion of a large sparse matrix, subject to non-negativity constraint. As a challenging prototype of this problem, we have in mind the well-known Netflix problem. Our method is based on the derivation of a constrained gradient system and its numerical integration. The methods we propose are based on the constrained minimization of a functional associated to the low-rank completion matrix having minimal distance to the given matrix. In the main 2-level method, the low-rank matrix is expressed in the form of the non-negative factorization X = eUVT, where the factors U and V are assumed to be normalized with unit Frobenius norm. In the inner level—for a given e—we minimize the functional; in the outer level, we tune the parameter e until we reach a solution. Numerical experiments on well-known large test matrices show the effectiveness of the method when compared with other algorithms available in the literature.

Rank deficient Hankel matrices are at the core of several applications. However, in practice, the coefficients of these matrices are noisy due to e.g. measurements errors and computational errors, so generically the involved matrices are full rank. This motivates the problem of Hankel structured low-rank approximation. Structured low-rank approximation problems, in general, do not have a global and efficient solution technique. In this paper we propose a local optimization approach based on a two-levels iteration. Experimental results show that the proposed algorithm usually achieves good accuracy and shows a higher robustness with respect to the initial approximation, compared to alternative approaches.

For the computation of the lower spectral radius of a finite family of matrices that shares an invariant cone, two recent papers by Guglielmi and Protasov [8] and Guglielmi and Zennaro [9] make use of so-called antinorms. Antinorms are continuous, nonnegative, positively homogeneous and superadditive functions defined on the cone and turn out to be related to the lower spectral radius of the family in a similar way as norms are related to the joint spectral radius. In this paper, we revisit the theory of antinorms in a systematic way, filling in some theoretical holes, correcting a common mistake present in the literature and adding some new properties and results. In particular, we prove that, under suitable assumptions, the lower spectral radius is characterized by a Gelfand type limit computed on an antinorm.

International audience; We study asymptotic stability of continuous-time systems with mode-dependent guaranteed dwell time. These systems are reformulated as special cases of a general class of mixed (discrete-continuous) linear switching systems on graphs, in which some modes correspond to discrete actions and some others correspond to continuous-time evolutions. Each discrete action has its own positive weight which accounts for its time-duration. We develop a theory of stability for the mixed systems; in particular, we prove the existence of an invariant Lyapunov norm for mixed systems on graphs and study its structure in various cases, including discrete-time systems for which discrete actions have inhomogeneous time durations. This allows us to adapt recent methods for the joint spectral radius computation (Gripenberg's algorithm and the Invariant Polytope Algorithm) to compute the Lyapunov exponent of mixed systems on graphs.

We consider the problem of computing (approximate) greatest common divisors for matrix polynomials and we present some related facts and applications in system and control theory. The main application is to compute the distance of a controllable multi-input multi-output system to the set of uncontrollable ones; then we describe some related results. Index Terms— Matrix polynomial common factor - Distance to uncontrollability - Behavioral approach

We consider linear dynamical systems with a structure of a multigraph. The vertices are associated to linear spaces and the edges correspond to linear maps between those spaces. We analyse the asymptotic growth of trajectories (associated to paths along the multigraph), the stability and the stabilizability problems. This generalizes the classical linear switching systems and their recent extensions to Markovian systems, to systems generated by regular languages, etc. We show that an arbitrary system can be factorized into several irreducible systems on strongly connected multigraphs. For the latter systems, we prove the existence of invariant (Barabanov) multinorm and derive a method for its construction. The method works for a vast majority of systems and finds the joint spectral radius (Lyapunov exponent). Numerical examples are presented and applications to the study of fractals, attractors, and multistep methods for ODEs are discussed.

Piece-wise smooth differential equations (their regularization, numerical integration, and classification of solutions) is the topic of the present work. The behaviour close to one discontinuity surface and also the entering into the intersection of two discontinuity surfaces is well understood. Here, we study the solutions that exit a codimension-2 sliding mode. Some results are expected, others come as a surprise. We are able to explain situations, where difficulties in numerical computations are reported in the recent literature. The analysis is based on asymptotic expansions for singularly perturbed problems and on the study of a time-parameterized two-dimensional dynamical system (hidden dynamics). Various situations are illustrated by examples.

The problem of computing the distance of two real coprime polynomials to the set of polynomials with a nontrivial greatest common divisor (GCD) appears in computer algebra, signal processing, and control theory. It has been studied in the literature under the names approximate common divisor, $\varepsilon$-GCD, and distance to uncontrollability. Existing solution methods use different types of local optimization methods and require a user-defined initial approximation. In this paper, we propose a new method that allows us to include constraints on the coefficients of the polynomials. Moreover, the method proposed in the paper is more robust to the initial approximation than the Newton-type optimization methods available in the literature. Our approach consists of two steps: (1) reformulate the problem as a problem of determining the structured distance to singularity of an associated Sylvester matrix, and (2) integrate a system of ODEs, which describes the gradient associated to the functional to be minim...

We consider the problem of stabilizing a matrix by a correction of minimal norm: Given a square matrix that has some eigenvalues with positive real part, find the nearest matrix having no eigenvalu...

A novel method for approximating structured singular values (also known as $\mu$-values) is proposed and investigated. These quantities constitute an important tool in the stability analysis of uncertain linear control systems as well as in structured eigenvalue perturbation theory. Our approach consists of an inner-outer iteration. In the outer iteration, a Newton method is used to adjust the perturbation level. The inner iteration solves a gradient system associated with an optimization problem on the manifold induced by the structure. Numerical results and comparison with the well-known MATLAB function mussv, implemented in the MATLAB Control Toolbox, illustrate the behavior of the method.

Computation of (approximate) polynomials common factors is an important problem in several fields of science, like control theory and signal processing. While the problem has been widely studied for scalar polynomials, the scientific literature in the framework of matrix polynomials seems to be limited to the problem of exact greatest common divisor computation. In this paper, we generalize two algorithms from scalar to matrix polynomials. The first one is fast and simple. The second one is more accurate but computationally more expensive. We test the performances of the two algorithms and observe similar behavior to the one in the scalar case. Finally we describe an application to multi-input multi-output linear time-invariant dynamical systems.

Computing the greatest common divisor of a set of polynomials is a problem which plays an important role in different fields, such as linear system, control, and network theory. In practice, the polynomials are obtained through measurements and computations, so that their coefficients are inexact. This poses the problem of computing an approximate common factor. We propose an improvement and a generalization of the method recently proposed in Guglielmi et al. (SIAM J. Numer. Anal. 55, 1456–1482, 2017), which restates the problem as a (structured) distance to singularity of the Sylvester matrix. We generalize the algorithm in order to work with more than 2 polynomials and to compute an Approximate GCD (Greatest Common Divisor) of degree k ≥ 1; moreover, we show that the algorithm becomes faster by replacing the eigenvalues by the singular values.

As an indicator of the stability of spectral clustering of an undirected weighted graph into $k$ clusters, the $k$th spectral gap of the graph Laplacian is often considered. The $k$th spectral gap is characterized in this paper as an unstructured distance to ambiguity, namely as the minimal distance of the Laplacian to arbitrary symmetric matrices with vanishing $k$th spectral gap. As a conceptually more appropriate measure of stability, the structured distance to ambiguity of the $k$-clustering is introduced as the minimal distance of the Laplacian to Laplacians of graphs with the same vertices and edges but with weights that are perturbed such that the $k$th spectral gap vanishes. To compute a solution to this matrix nearness problem, a two-level iterative algorithm is proposed that uses a constrained gradient system of matrix differential equations in the inner iteration and a one-dimensional optimization of the perturbation size in the outer iteration. The structured and unstructured distances to ambiguity are compared on some example graphs. The numerical experiments show, in particular, that selecting the number $k$ of clusters according to the criterion of maximal stability can lead to different results for the structured and unstructured stability indicators.

Under the direction of the Institute for Informatics and Telematics in Pisa, Calcolo publishes original contributions on numerical analysis and its applications, and on the theory of computation. The main focus of the journal is on numerical linear algebra, approximation theory and its applications, numerical solutions of differential and integral equations, computational complexity, algorithmics, mathematical aspects of computer science, and optimization theory. Calcolo contains expository papers which introduce emerging topics as well as a abstracts of PhD theses, news and reports from conferences, and book reviews.

A linear gyroscopic system is of the form: $$\displaystyle M \ddot x + G\dot x + K x = 0, $$ where the mass matrix M is a symmetric positive definite real matrix, the gyroscopic matrix G is real and skew symmetric, and the stiffness matrix K is real and symmetric. The system is stable if and only if the quadratic eigenvalue problem \(\det (\lambda ^2 M+\lambda G + K)=0\) has all eigenvalues on the imaginary axis.

In this paper, we study scalar multivariate non-stationary subdivision schemes with integer dilation matrix M and present a unifying, general approach for checking their convergence and for determining their Hölder regularity (latter in the case \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M = mI, m \ge 2$$\end{document}M=mI,m≥2). The combination of the concepts of asymptotic similarity and approximate sum rules allows us to link stationary and non-stationary settings and to employ recent advances in methods for exact computation of the joint spectral radius. As an application, we prove a recent conjecture by Dyn et al. on the Hölder regularity of the generalized Daubechies wavelets. We illustrate our results with several examples.

We propose a method for computing the distance of a stable polynomial to the set of unstable ones (both in the Hurwitz and in the Schur case). The method is based on the reformulation of the problem as the structured distance to instability of a companion matrix associated to a polynomial. We first introduce the structured e -pseudospectrum of a companion matrix and write a system of ordinary differential equations which maximize the real part (or the absolute value) of elements of the structured e -pseudospectrum and then exploit the knowledge of the derivative of the maximizers with respect to e to devise a quadratically convergent iteration. Furthermore we use a variant of the same ODEs to compute the boundary of structured pseudospectra and compare them to unstructured ones. An extension to constrained perturbations is also considered.

We generalize the recent invariant polytope algorithm for computing the joint spectral radius and extend it to a wider class of matrix sets. This, in particular, makes the algorithm applicable to sets of matrices that have finitely many spectrum maximizing products. A criterion of convergence of the algorithm is proved. As an application we solve two challenging computational open problems. First we find the regularity of the Butterfly subdivision scheme for various parameters $\omega$. In the “most regular” case $\omega = \frac{1}{16}$, we prove that the limit function has Holder exponent 2 and its derivative is “almost Lipschitz” with logarithmic factor 2. Second we compute the Holder exponent of Daubechies wavelets of high order.

In a recent paper, Rostami [SIAM J. Sci. Comput, 37 (2015), pp. S447--S471] has presented an interesting algorithm for the computation of the real pseudospectral abscissa and the real stability radius (aka the distance to instability) of a square matrix $A \in \mathbb R^{n,n}$ in the spectral norm. The method is particularly well suited for large sparse matrices. The algorithm to compute the stability radius is based on the alternation of an iterative method in the rank-$2$ manifold of real matrices and a Lyapunov inverse iteration. The algorithm shows quadratic convergence in all experiments, but this peculiar property is only conjectured. In this paper, we provide a rigorous proof and propose a modification which replaces the Lyapunov inverse iteration and speeds up the method. Moreover, we interpret the algorithm to compute the real $\varepsilon$-pseudospectral abscissa as a nonstandard discretization of the ODE-based method given in Guglielmi and Lubich [SIAM J. Matrix Anal. Appl., 34 (2013), pp. 40--...

We present a new approach for constructing polytope Lyapunov functions for continuous-time linear switching systems (LSS). This allows us to decide the stability of LSS and to compute the Lyapunov exponent with a good precision in relatively high dimensions. The same technique is also extended for stabilizability of positive systems by evaluating a polytope concave Lyapunov function ("antinorm") in the cone. The method is based on a suitable discretization of the underlying continuous system and provides both a lower and an upper bound for the Lyapunov exponent. The absolute error in the Lyapunov exponent computation is estimated from above and proved to be linear in the dwell time. The practical efficiency of the new method is demonstrated in several examples and in the list of numerical experiments with randomly generated matrices of dimensions up to $10$ (for general linear systems) and up to $100$ (for positive systems). The development of the method is based on several theoretical results proved in the paper: the existence of monotone invariant norms and antinorms for positively irreducible systems, the equivalence of all contractive norms for stable systems and the linear convergence theorem.

In this paper a novel contour integral method is proposed for linear convection-diffusion equations. The method is based on the inversion of the Laplace transform and makes use of a contour given by an elliptic arc joined symmetrically to two half-lines. The trapezoidal rule is the chosen integration method for the numerical inversion of the Laplace transform, due to its well-known fast convergence properties when applied to analytic functions. Error estimates are provided as well as careful indications about the choice of several involved parameters. The method selects the elliptic arc in the integration contour by an algorithmic strategy based on the computation of pseudospectral level sets of the discretized differential operator. In this sense the method is general and can be applied to any linear convection-diffusion equation without knowing any a priori information about its pseudospectral geometry. Numerical experiments performed on the Black-Scholes ($1D$) and Heston ($2D$) equations show that the method is competitive with other contour integral methods available in the literature.

We consider the problem of computing the greatest common divisor of a set of univariate polynomials and present applications of this problem in system theory and signal processing. One application is blind system identification: given the responses of a system to unknown inputs, find the system. Assuming that the unknown system is finite impulse response and at least two experiments are done with inputs that have finite support and their Z-transforms have no common factors, the impulse response of the system can be computed up to a scaling factor as the greatest common divisor of the Z-transforms of the outputs. Other applications of greatest common divisor problem in system theory and signal processing are finding the distance of a system to the set of uncontrollable systems and common dynamics estimation in a multi-channel sum-of-exponentials model.

We consider the problem of computing the closest stable/unstable non-negative matrix to a given real matrix. This problem is important in the study of linear dynamical systems, numerical methods, etc. The distance between matrices is measured in the Frobenius norm. The problem is addressed for two types of stability: the Schur stability (the matrix is stable if its spectral radius is smaller than one) and the Hurwitz stability (the matrix is stable if its spectral abscissa is negative). We show that the closest unstable matrix can always be explicitly found. For the closest stable matrix, we present an iterative algorithm which converges to a local minimum with a linear rate. It is shown that the total number of local minima can be exponential in the dimension. Numerical results and the complexity estimates are presented.

Ordinary differential equations with discontinuous right-hand side, where the discontinuity of the vector field arises on smooth surfaces of the phase space, are the topic of this work. The main emphasis is the study of solutions close to the intersection of two discontinuity surfaces. There, the so-called hidden dynamics describe the smooth transition from ingoing to outgoing solution directions, which occurs instantaneously in the jump discontinuity of the vector field. This paper presents a complete classification of such transitions (assuming the vector fields surrounding the intersection are transversal to it). Since the hidden dynamics are realized by standard space regularizations, much insight is obtained for them. One can predict, in the case of multiple solutions of the discontinuous problem, which solution (classical or sliding mode) will be approximated after entering the intersection of two discontinuity surfaces. A novel modification of space regularizations is presented that permits us to a...

We propose a fast method to approximate the real stability radius of a linear dynamical system with output feedback, where the perturbations are restricted to be real valued and bounded with respect to the Frobenius norm. Our work builds on a number of scalable algorithms that have been proposed in recent years, ranging from methods that approximate the complex or real pseudospectral abscissa and radius of large sparse matrices (and generalizations of these methods for pseudospectra to spectral value sets) to algorithms for approximating the complex stability radius (the reciprocal of the $H_\infty$ norm). Although our algorithm is guaranteed to find only upper bounds to the real stability radius, it seems quite effective in practice. As far as we know, this is the first algorithm that addresses the Frobenius-norm version of this problem. Because the cost mainly consists of computing the eigenvalue with maximal real part for continuous-time systems (or modulus for discrete-time systems) of a sequence of matrices, our algorithm remains very efficient for large-scale systems provided that the system matrices are sparse.

Given a regular matrix pencil $A + \mu E$, we consider the problem of determining the nearest singular matrix pencil with respect to the Frobenius norm. We present new approaches based on the solut...

We present several results describing the interplay between the max algebraic joint spectral radius (JSR) for compact sets of matrices and suitably defined matrix norms. In particular, we extend a classical result for the conventional algebra, showing that the max algebraic JSR can be described in terms of induced norms of the matrices in the set. We also show that for a set generating an irreducible semigroup (in a cone-theoretic sense), a monotone Barabanov norm always exists. This fact is then used to show that the max algebraic JSR is locally Lipschitz continuous on the space of compact irreducible sets of matrices with respect to the Hausdorff distance. We then prove that the max algebraic JSR is locally Hoelder continuous on the space of compact sets of nonnegative matrices. Finally, we prove a strict monotonicity property for the max algebraic JSR that echoes a fact for the classical JSR. The single matrix characterisation of the max algebraic JSR plays a vital role in our proofs.

We study differential equations that lead to extremal points in symplectic pseudospectra. In a two-level approach, where on the inner level we compute extremizers of the symplectic $\varepsilon$-pseudospectrum for a given $\varepsilon$ and on the outer level we optimize over $\varepsilon$, this is used to solve symplectic matrix nearness problems such as the following: For a symplectic matrix with eigenvalues of unit modulus, we aim to determine the nearest complex symplectic matrix such that some or all eigenvalues leave the complex unit circle. Conversely, for a symplectic matrix with all eigenvalues lying off the unit circle, we consider the problem of computing the nearest symplectic matrix that has an eigenvalue on the unit circle.

This article presents two regularization techniques for systems of state dependent neutral delay differential equations which have a discontinuity in the derivative of the solution at the initial point. Such problems have a rich dynamics and besides classical solutions can have weak solutions in the sense of Utkin. Both of the presented techniques permit the numerical solution of such problems with the code {\sc RADAR5} which is designed to compute classical solutions of stiff and differential algebraic (state dependent) delay equations.

For differential equations with discontinuous right-hand side and, in particular, for neutral delay equations it may happen that classical solutions do no exist beyond a certain time instant. In this situation, it is common to consider weak solutions of Utkin (Filippov) type. This article extends the concept of weak solutions and proposes a new regularization which eliminates the discontinuities. Codimension- 1 and codimension- 2 weak solutions are considered. Numerical experiments show the advantages of the new regularization.

This talk discusses the concept of solutions of piecewise smooth, discontinuous differential equations. Recent results of [1] are presented, which give new insight into the solution close to the intersection of two discontinuity hyper-surfaces. The study is connected to space regularizations where, close to the discontinuities, the vector field is smoothed out by linear (bi-linear) interpolation.