Question: Checking the validity of output of dsolve, numeric

Dear mapleprimes ODE enthusiasts,

I am analyzing a stiff 3D system characterizes by fast-slow dynamics. I believe I can show analytically that, for some parameter values, in a neighborhood of some critical points (where the system is stationary), the dynamic system has a local center manifold that is center-stable. I have simulated the system for different parameter values.I am interested in one particular trajectory. With a random set of initial values, the system will typically collapse onto itself or explode, but for specific initial values (to be determined) I believe it is possible to follow a trajectory that converges to the center-stable manifold and thus to the given critical point (proof of the center-stable manifold is local, but I'm hoping/assuming that it might be global). To sum up, I have a two-point boundary value problem, where I fix one end-point (the critical point, that which is stationary) but leave the other end-point free.

I have found it difficult to find such a center-stable trajectory for some parameter values. I proceeded in two steps: Step 1: dsolve the center-stable manifold approximation (it can actually be solved exactly, so my use of dsolve,numeric in that part of the worksheet could be avoided but it's typically fast so I kept it), Step 2: dsolve the non-approximated boundary value problem with end-point conditions near the critical point specified to coincide with the center manifold approximation.

I have managed to simulate the center-stable trajectory for some parameter values, but have run into difficulty with others. In some cases, I have failed to get convergence. In other cases, I have succeeded in getting convergence, but the resulting trajectory looks nothing like the center manifold approximation and is therefore suspect. (However, since the center manifold approximation is local, a visual discrepancy away from a neighborhood of the critical point is not impossible in theory) Hence my question: Do you have any suggestions on how I could check that my simulation is approximately correct.

Any other comments on other aspects of the code are always welcome.

Because the system and the calculations are long, I will not attempt to copy-paste the code here. Instead, I attach a worksheet that summarizes some of my calculations using dsolve. I simulate the system with two sets of parameter values (named parameters[1] and parameters[10], to make the worksheet legible I have deleted the other 8 sets of values in between as well as other intermediate calculations).

To sum up what I do: I define the dynamic system (x,c,q), I then define a transformed system (u,r,q), I then define a center manifold approximation of the transformed system. I index the critical point I'm interested in by the stationary value of q, denoted qss. From this value follow all the other stationary values of interest, i.e. xss=ifx(qss), css=f(xss), rss=qss, and uss=0, where ifx and f are functions defined early on in the worksheet.

I plot the resulting solutions in the (x,c) plane. The orange line is the trajectory that converges to the critical point (xss,css,qss). The dotted red trajectory is the center manifold approximation. They should coincide near a neighborhood of the critical point.

I have found it typically easier to simulate a trajectory from below (x<xss) than from above (x>xss). I have experimented with high levels of digits and error intolerance. I have set up events to halt integration when the variables go off regions of interest/stability (to avoid unnecessary computations). I have experimented with time running forwards and backwards. Here's a typical call to dsolve:

### Solution OFF Center Manifold, Running Time FORWARD
Digits := 100:
Precision := 12:                          ## Tweak This. High precision needed if r(t) close to 0.
AbsErr := 10.^(-Precision):
RelErr := 10.^(-Precision):
T := 100:                               ## T is Positive
INI := op([ u(T) = uT, r(T) = rT, q(T) = qT ]): 'INI' = evalf(INI);  ## Start-Point Conditions
SOL := dsolve(
  [SYS, INI]
  , VAR
  , range = 0 .. T
  , type = numeric
  , abserr = AbsErr
  , relerr = RelErr
  , output = listprocedure
  , stiff = true
  , events = [ [r(t) = -0.1 .. qss, halt], [u(t) = -css .. 0, halt] ] ## Tweak This
);

With parameters[1] I obtain this trajectory:

The trajectory is consistent with what is known about these kind of trajectories in related models (but this particular model being different in some ways, this is little comfort). The slope of the trajectory (orange line) looks a little off the slope of the center manifold approximation (red dots). But this could be because the approximation is only valid in a small neighborhood.

with parameters[10] I obtain this trajectory:

Here, what worries me is the simulated trajectory (orange line) above the critical value xss. It has a hump-shape and looks quite different from the center manifold approximation (red dots). Can it be right? Is there a way to check?

Mapleprimes-CenterMa.mwMapleprimes-CenterMa.mw