Numerical integration "backwards"

Hello everyone,

I have a system of nonlinear equations on the form

x1' = f1(x1,x2,x3)

x2' = f2(x1,x2,x3)

x3' = f3(x1,x2,x3)

Let the flow F(x,t) be the solution to the system, that is, F(x0,t) as a function of t is an integral curve that goes through x0 at t=0.

I need the value of the inverse flow to some particular values of x0 and t0, that is, the value of F^-1(x0,t0) = F(x0,-t0).

Wen I use Matlab's odeXX functions to integrate F() numerically, I get well-behaved integral curves. Now, the problem is, when I try to use them to integrate F^-1() by reversing the time span parameter (or by exchanging the signs of the right hand sides of the differential equations), then apparently a numerical unstability occurs, and the integral curve grows extremely fast.

Does any one know how I can get the curve for F^-1(x0,t0), for particular t0 and x0?

Code used so far:

%code to compute the solution for F(x0,t0):
[t1,flow] = ode15s(@f1f2f3,[0 t0],x0);  

%code to compute the solution for F^-1(x0,t0)=F(x0,-t0):
[t2,flowinv] = ode15s(@f1f2f3,[t0 0],x0);

[Teja Muppirala]

It'd be interesting to see the actual equations though.

Answer 1

[Andrew Newell]

You seem to be talking about equations for some physical phenomenon like fluid flow. In general, you can't integrate physical equations backwards. You can only do it if you don't have dissipative terms like friction. The numerical blowup is probably telling you that you have such a term. This is not really a MATLAB issue.

EDIT: You could get an intuitive feel for whether it is feasible or not by plotting a phase portrait of the system, i.e., choose a grid of reasonable starting values, calculate the flow, and then plot all the curves together. Do you have attractors?

[Teja Muppirala]

I wouldn't go so far as to say that it's ALWAYS impossible to integrate backwards when there are dissipative terms. For example, here is a damped pendulum, integrated both forwards and backwards.

mydx = @(t,x) [x(2); -0.1*x(2)-x(1)];

[t_f,X_f] = ode45(mydx,[0:.1:10],[1; 1])

[t_b,X_b] = ode45(mydx,[0:-.1:-10],[1; 1])

plot(t_f,X_f)

hold all

plot(t_b,X_b)

[José Goulart]

Thank you for your comments.

Indeed, I know that this situation is caused by some properties of my particular problem itself. However, I was wondering if there is a specific routine/function to handle situations like this, using some sort of method which prevents this blowup to occur.

But, according to what you wrote, I suppose that unfortunately the answer is no.

[José Goulart]

Yes, Andrew, I have only one attractor, which is the fixed point at the origin.

[Andrew Newell]

Then you might be able to run the simulation backwards as long as you are far enough from the origin.

[bym]

@Andrew

Is a phase portrait different from a Poincaré map?

[Andrew Newell]

Yes. The phase portrait is just a picture of flow lines while the Poincaré map is a less intuitive plot that looks at intersections of trajectories in some traversal section in phase space. Sorry, I can't really explain it in a comment. The point is that if the orbits are periodic they'll keep intersecting at the same points, but if it's chaotic you'll get a much more complex pattern.

[José Goulart]

Even when I start the simulation in a state x0 very far from the origin (picking state components with even rather unrealistic values for the quantities they represent), the computation of the time-reversed trajectory only works correctly for a very small interval of time (some miliseconds). After this, a numerical blowup still occurs.

Thinking about your previous comment, I'm starting to think it will be really difficult to get this time-reversed trajectory for my system, because the equations have very strong dissipative terms. Computing the conventional (forward) trajectory numerically, I see that even for initial states quite far from the origin, the "decaying time" (i.e., the time required for the trajectory to get very close to the origin) is quite short (the observed times were never longer than about 30 ms).

Moreover, when I compute the Jacobian of the equations around the initial state, its eigenvalues have a really large magnitude. So, I conclude that by inverting the signs of the equations, I have a very unstable set of equations that would need and infinitesimal step size to be numerically integrated.