Using ode45 with terminal values

[Alexandra Loubeau]

I am trying to use ode45 to solve an IVP problem with terminal values
(for example to step backwards from t=15 to t=-15). I have seen
several posts on the MATLAB Newsgroup claiming that this is possible
(by declaring a time interval [15 -15], but I haven't gotten it to
work yet. Any help would be greatly appreciated!

Here is an example where I am recreating a function from MATLAB Help
(Example: Solving an IVP ODE (van der Pol Equation, Nonstiff) ). I
can solve it correctly in the forward direction, but the result is
incorrect for backwards integration.

%Forward direction works
[t,y] = ode45(@vdp1,[0 20],[2; 0]);

figure;
plot(t,y(:,1),'-',t,y(:,2),'--')
title('Solution of van der Pol Equation, \mu = 1');
xlabel('time t');
ylabel('solution y');
legend('y_1','y_2','Location','Best')

%Backwards integration doesn't work
%periodic function, so use same values for terminal %conditions
[t,y] = ode45(@vdp1,[20 0],[2; 0]);

figure;
plot(t,y(:,1),'-',t,y(:,2),'--')
title('Solution of van der Pol Equation, \mu = 1');
xlabel('time t');
ylabel('solution y');
legend('y_1','y_2','Location','Best')

%van der Pol equation
function dydt = vdp1(t,y)
dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];

[helper]

If you think the backward integration is incorrect, rerun your
forward integration using the final conditions of the backward
integration "y(end,2)" as the initial conditions.

You will see the backward and forward integration is identical. I'm
sure someone else can better explain why this is happening.

[Alexandra Loubeau]

Thanks for your reply. Yes, if I use the y(end,:) values from the
backward integration as the inital conditions for the forward
integration, I get the same results. But why doesn't this work the
other way? If you take the y(end,:) values from the forward
integration and use them as the terminal conditions for the backward
integration, I don't get the same answer. By the way, the original
forward integration with initial conditions [2;0] yields the same
plot as in the MATLAB Help (while the backward integration does not
match). Thanks.

[Jacek Kierzenka]

Alexandra,

Let me help sorting this out:

1) With MATLAB ODE solvers, numerical integration is possible in both
directions.
Take something simple, like y' = 1, y(0) = 0, y(1)=1.
Using anonymous functions, this could be coded as
>> [t0,y0] = ode45(@(t,y) 1, [0,1],0);
>> [t1,y1] = ode45(@(t,y) 1,[1,0],1);
>> plot(t0,y0,t1,y1,'o')
You will see that the plots overlap.

2) Due to problem stability, integrating ODEs in one direction could
sometimes be _much_ more difficult than integrating them in the opposite
direction. Basically, numerical errors are damped when integrating in one
direction, but are amplified when integrating in the other. Apparently,
VDP1 in 'positive' direction is stable, but in 'negative' direction is not.
This is why you were able to start your 'positive' integration with the
'negative' solution and match the 'negative' results, but not the other way
around.

All that said, if you do your computations asking for enough accuracy, the
results will match (well, to the naked eye, at least). Try this:

>> options = odeset('RelTol',1e-9,'AbsTol',1e-12);
>> [t0,y0] = ode45(@vdp1,[0, 20],[2;0],options);
>> [t1,y1] = ode45(@vdp1,[20,0],y0(end,:),options);
>> plot(t0,y0,t1(1:20:end),y1(1:20:end,:),'o')

You will see that the 'negative' solution starts diverging close to t=0.
Asking for even more accuracy will get it back on track...

Hope that explains it.
-Jacek

-------------------------------------

"Alexandra Loubeau" <alou...@psu.NO.SPAM.edu> wrote in message
news:ef27...@webx.raydaftYaTP...

[Alexandra Loubeau]

Thanks, Jacek. Tightening error tolerances seems to do the trick
with the van der Pol example. Now I just have to apply this method
to my actual problem...

Thanks!

Alexandra

[helper]

helper wrote:

> I'm sure someone else can better explain why this is happening.

See, I was right! Thanks Jacek.