VIRTUAL SLOW MANIFOLDS

1. VIRTUAL SLOW MANIFOLDS C. W. Gear with I. Kevrekidis & D. Givon Princeton

2. What is a virtual slow manifold? • Some problems exhibit continuing fast behavior but • the fast behavior of interest can be characterized by a few slow variables, and • the local fast behavior can be generated very quickly from the differential system from one single starting point (there will typically be many possible starting points that will generate essential the same fast behavior). • Examples: • An oscillatory problem with a slowly changing wave shape (and "period") • Some problems with multiple very fast oscillations • Some problems that are stochastic in some of the variables A virtual slow manifold is a manifold each of whose points is a one-one mapping of each local fast solution of a problem with the above characteristics.

3. For the manifold to be a useful construct, it must be easy to compute this mapping in either direction, i.e., to find an approximation to the point on the virtual manifold from a "short" computation, and to compute an initial value that will give a solution corresponding to a point on the manifold. • We also need to be able to compute the derivative – or chord slope – of the solution on the slow manifold. This can be done with the following steps: • Given a point on the virtual slow manifold compute an initial value of the full system (lifting) • Integrate for some time and from the solution compute one or more additional points on the slow manifold (restriction) • Compute the chord slope from two of these points (or other descriptors) • This process can then be used in conventional numerical algorithms to find "steady states" with a zero finder, to integrate forward in time over long periods using projective integration, compute bifurcation diagrams, etc.

4. Plan of talk • Discuss a trivial oscillatory example • Consider simple sample multiple oscillations and stochastic systems • Look at the stochastic problem in some algorithmic detail to understand some of the issues • Look at the multiple oscillator problem to examine additional issues

5. Example: Nearly periodic problem A slow manifold on this is simple – a smooth curve that passes through points about one period apart. If the problem is autonomous, a point on this manifoldisan initial value for a solution

6. These were called "quasi envelopes" in Petzold's 1978 thesis. Part of what we are doing could be viewed as a generalization of this.

7. Integrate one cycle calculate chord slope and execute projective step The virtual slow manifold is a set of points that represent all possible "local" solutions, so the evolution of a fast solution is a path in the virtual slow manifold. Our objective is, of course, to find a way to compute this path without computing the solution over the fast problem for the whole interval of interest. In Petzold's thesis the chord over approximately a period was calculated and used in an integration formula – using what we now call projective integration.

8. Integrate two cycles calculate chord slope from last cycle and execute projective step Recently we have been looking at some problems that also include "stiff" components so we integrate through one (or more) cycles before computing the chord slope to damp the stiff components.

9. The selected points do not have to be precisely periodic to give a "slow" manifold – magenta line

10. And the virtual slow manifold is not unique – any of the three shown satisfy our definition but all of these virtual slow manifolds have the property that each point on them generates a solution (assuming the system is autonomous). In practical computation we adopt some definition that makes the manifold unique.

11. If and are not large - e.g., O(1) – and we are not interested in phase information, then this has a 3-dimensional "virtual slow manifold" characterized by z and the amplitudes of the two oscillators (z has a a small oscillatory behavior of order ). For each that is large, the amplitude of the corresponding oscillator approaches z very quickly, and the manifold is of one dimension less. Now consider a simple multiple oscillator example:

12. Simple Stochastic DE Consider a fast-slow system that is stochastic in the fast dimensions: • Here, ε has to be essentially zero so that • the stochastic part reaches it's locally steady state distribution very rapidly and • the slow components do not exhibit significant stochastic behavior over the interval of interest. • We assume that our problem has this fast-slow behavior, but that our observable variables, w, are unknown transformations so that we don't have direct access to the fast and slow components.

13. What problem description are we working with? • We almost never have the equations, only a code that will evaluate something • The legacy code issue (although this is less common) • The system is large and described in a modeling language – for example, as an electrical network, an integrated circuit layout description, a Simulink description, … • While in there is, in principle, a set of equations that might be obtained, that step is usually not practical. • Usually all we can obtain from the code is the solution at a future point (or set of points) at some time spacing (h) that we may be able to specify. • In the case of stochastic equations, the inner code mechanism is often not an SDE but a random number generator (the SDE is seen as convenient mathematical description – but it is often not a good simulation tool)

14. Stochastic problems If we obtained the output at a set of time points separated by an extremely small h (but large compared to ε) the computed solution would be a "cloud" of points almost in a low-dimensional manifold approximating a probability density function. If we ran with a somewhat larger h the cloud would "thicken" but still provide approximate information about the probability density of the solution locally and from this it should be possible to extract some information about the slow behavior of the system. At first one might think of trying to compute the average derivative over the cloud (a method used for fast-slow systems by Vanden-Eijnden) but is does not work when we do not know the slow components. Instead, we have to look for characteristic(s) that is(are) slowly changing, for example, the cloud average. Clouds at two well separated times.

15. Example: This example has p = 2 and q = 1; i.e., two ODEs coupled to one fast stochastic system and then subject to a non-linear transformation. In this example, we show the trajectory of the average for a particular solution and the clouds at isolated places. Here we could define a virtual slow manifold from all locally averaged solutions (or similar statistic).

16. MS We want to characterize clouds by some computable value(s), P(cloud), such as the average – the black dot in the figure. If we use the average, P, it exists in the same space as the problem but may not lie in the cloud. If we were to perform this operation for all possible starting values we would get a collection of P's which would lie on a p-dimensional manifold, MS (2-dimensional in this case).

17. Our objective is to use the same type of approach as in the simple oscillatory problem: Run for a while to compute two points on the virtual slow manifold; use extrapolation (illustrated here as first order) compute the chord; But, in the extrapolation process we are estimating a point on the virtual slow manifold, whereas we need to find the corresponding cloud to continue the process.

18. Estimated average point Integration starting point Orthogonal projection onto slow manifold There is another issue to deal with – the projected average point typically does not lie on the p-dimensional virtual slow manifold so we would like to know the local tangent plane to the slow manifold and find an integration starting point that gets us to the orthogonal projection of the estimated average point.

19. Numerical integration to get a "cloud" of points (approximating the PDF) will advance the solution in time, so the cloud we get will not represent the instantaneous PDF of the solution but an averaged value of the PDFs over the time interval used for averaging. However, this does not matter. There still exists a q-dimensional manifold (a "spine") of starting values that will give (asymptotically) the same cloud, and all we need to do is to find any starting point on this manifold which leads to a cloud average that matches the desired average in the (approximate)direction of the slow manifold.

20. Computing a "cloud" • Select a starting point. Simulate over a few intervals (to remove any "stiff" components) then run for N more intervals to compute the cloud statistic (average) • Computing a numerical solution: • Compute the first cloud. • Continue simulation to compute a second cloud • Restrict to the solution average and compute the chord • Use projective integration in the slow manifold. • Assuming(A) we have an approximation to the tangent plane of the virtual slow manifold, Liftto find(B) a starting point for a cloud integration that yields a statistic (average) that has the same projection in the tangent plane. • Update(C) the tangent plane approximation – using any information gathered during steps 2 and 4. • We have computed a cloud at the end of step 4, so we can use this as the first cloud in the next step and return to step 2. • How to handle (A), (B), and (C)?

21. (A) Finding an approximation to the tangent plane to the slow manifold Initially it seems that we will have to do some "exploration" of the solution numerically (unless we have some prior knowledge). Presumably we know the total dimension n = p + q but we may not know p the dimension of the slow manifold. The direct way to find p is to compute clouds corresponding to n perturbations of the starting point (but this could be very expensive). If we can find the dimension, q, of the first cloud it would be less expensive because we then know p. A possible approach here is the use of diffusion maps (R. R. Coifman et al. "Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps," PNAS 2005 102: 7426-7431) We might also be able to get an approximation to the tangent manifold to the cloud, MC, at a specific cloud point using a diffusion map (this could be valuable in task (B) -- solving for the starting point of a cloud integration) but it will be particular valuable for estimating the slow manifold tangent plane, MS, because we can limit ourselves to p perturbations orthogonal to MC. (We’ll discuss this later)

22. Point intersection Tangent direction planes Non-point intersection (B) Finding a starting point for a cloud integration We are given the slow manifold tangent direction and the extrapolated value of the average. We want to find a starting point for a cloud integration such that the projections of the cloud average and the extrapolated value agree. This is a p-dimensional non-linear equation which could in principle be solved with any standard package. However there is a difficulty that seems make this the most challenging part of the problem. The unknowns exists in an n = p + q dimensional space, so we would really like to restrict our search to a p-dimensional linear manifold – but we have to have one that only has a single point intersection with the spine of the desired cloud. It is tempting to consider a linear manifold that is parallel to the slow manifold tangent direction we have calculated but there is no guarantee that the spine of the cloud does not have more than a point intersection with such a manifold in some places

23. Starting point for iteration We need an initial point to begin an iterative process for solving for the integration starting value. The obvious choice is to extrapolate a value from prior starting points using a formula similar to the projection integration formula. We can then use a Newton-like method provided we have an approximation to the Jacobian relating perturbations in the integration starting value to perturbations in the average (both of these being in directions in the slow manifold tangent direction). Depending on the size of p we could opt to calculate an approximation to the Jacobian numerically from time to time (when convergence is slow) try to update it using Broyden-like corrections, or opt for matrix-free approaches.

24. When we get close to a point where the spine has a non-trivial intersection with the tangent plane we will run into convergence problems because the Jacobian will be getting close to singular. It might be appropriate to compute a new Jacobian at this point because this will also give us a better approximation to the tangent plane (and if the difficulty arose because of a poor tangent plane approximation, this will solve it). However, the reason we may be having trouble could be that we have moved to a point on the spine that is out in the tail (where the non-linearity of the transformation may give it a very different direction than near the "center"). therefore, we need to reset our iteration plane from time to time so that it is close to the “center” of the action. (We’ll discuss this later)

25. Rotated tangent direction MS (C) Updating the tangent plane approximation In the initial steps of projective integration we compute at least two points in the virtual slow manifold, so this gives a vector direction that should be approximately in the tangent plane. As we try to solve for the starting point of a cloud we will often iterate one or more times which will also give us vector directions approximately in the tangent plane. We can use these vectors to update the tangent direction data making as small a change as possible consistent with the new data New vector direction v Original tangent direction MS linear manifold in tangent plane orthogonal to v

26. Previous Example: (p = 2 and q = 1) Here we have computed two clouds (red and blues), used projective integration on their averages, and then solved for a starting point to find the next cloud. The two lines (black and red) are the projective solution and the solution obtained by simulation over the whole interval.

27. Finding the "center" of a cloud, estimating its dimension and finding an approximation to the tangent plane to it's spine. Finding a "center" is similar to a Google ranking algorithm. Define the "centerness" value of a point to be the weighted sum of the centerness values of nearby neighbors, where the weights decrease rapidly with distance – for example, where d is the distance and c is a normalizing scalar – actually the largest eigenvalue of the weight matrix. The "centerness" value of each point is then the value of the corresponding component of the eigenvector corresponding to the dominant eigenvalue. Example (next slide): 2000 points uniformly distributed in 2x2 x-y square, then assigned a z-value on a sphere of radius 1.2 and then perturbed in all three dimensions with Gaussian noise, standard deviation 0.3:

29. Estimating cloud dimension and a tangent plane: the diffusion map method Start with the weight matrix used before, but convert it to the transpose of a Markov matrix (a column stochastic matrix) by diving each column by its column sum. Suppose this matrix is M. Suppose that the coordinates of the N points in the cloud are represented as a d x N matrix A. The set of points represented by AM is the set of points we would get by replacing each point by a weighted sum of its nearby neighbors. This tends to move points near the outside of the cloud towards the center, but points in the interior that are fairly uniformly surrounded by neighbors on all sides are not moved significantly. Points on a curved lower-dimensional manifold tend to get "pulled into" the curve so that the curve is flattened out. If the operation is applied repeatedly to get AMmwe expect to see the cloud approach a lower-dimensional manifold – that is, its rank will reduce. We can possibly find the information only by studying the eigenstructure of M but the issues are not straightforward and it is hard to see how to get a robust algorithm of the sort we would want to put into a DE code. If we look at the SVD of the points after each transformation (translated so that their mean is zero) we seem to find a good separation in singular values. The next slide shows the result of five steps of the process (where the result has been shifted each time to keep the center point fixed):

30. 13.9478 13.5006 2.3726 29.3318 28.8486 15.6093 20.0675 19.6037 5.8131 4.7244 4.4271 0.1932 9.7205 9.3121 1.0141 6.7768 6.4218 0.4413

31. Singular Values of data

32. 1-D Case: If and are very large we can use z with any reasonable values of the other four variables as the one-dimensional virtual slow manifold. However, usually we will be dealing with a set of variables w = w(x1,x2,y1,y2,z) where the (non-linear) transformations are unknown, so we won't be able to able to find z directly from the solution. Multiple, fast, incommensurate oscillations situation:

33. This is the solution over about 200 cycles of the slower frequency:

34. While, in principle, one might be able to find a slow manifold that passes through points on solutions, it does not seem computationally feasible so we have to consider manifolds like a moving average – the red line below. The solution appears to "fill in" a region. In the plot below, we have shown a similar solution for occasional time bursts of about 200 cycles of the slowest frequency over a total time interval of the order of 106 cycles. Unfortunately we cannot find a "quasi-envelope" for such a problem (unless is the rational p/q with small p and q in which case we can use the smallest common integer multiple as the period).

35. As in the stochastic case, we have to consider functionals of the solution, for example, the "locally averaged" value of at least one of the variables in the 1-D case will contain enough information to determine the solution (to within a frequency dependent accuracy) but we must we must solve for an initial value as in the stochastic case. The local average will not tell us how to reconstruct the solution if we do not have the 1-D case because the averages maycontain no information about the amplitude of the oscillations. To get this we need to consider higher moments such as the local "variance" of the solution (e.g., a local average of the squared deviation from the local average). The virtual slow manifold now exists in a different space. For example, if we took the space of local averages and variances and small and we will find (approximately) a 3-dimensional slow manifold. We can use the same process: Restriction, low-dimensional numerical methods (such as projective integration) and lifting

36. (Preliminary) conclusions: It is possible to define low-dimensional manifolds for systems of the types discussed. If there is sufficient separation of time constants, computational savings are possible. The techniques are more likely to be attractive when we cannot realistically access the underlying equations.