ode45: can I reverse tspan for better initial conditions?

I'm using ode45 to solve/plot a 2nd-order differential equation in Matlab.

My tspan is from 0 to 0.25. But the initial condition near zero are ill-defined (Slope goes to infinity, complex values). Conditions at the end of tspan (near 0.25) are well defined (slope and value are zero).

Questions:

1. Can I reverse tspan, and use the "final conditions" as initial conditions?
2. Well, I know I can do it (see code below), and I get a plot that looks like what I expect, but is this a valid thing to do in general? Am I "getting lucky" in this one case?
3. ode45 provides numerical solutions and is not exact. Am I likely to have larger error after reversing tspan?

Here's my code (should run standalone):

    function ReverseTspan()
    % solve diff-eq backward from tspan end to tspan start using ode45()
    %  - Good initial conditions at the end, but not start.
    %  - Is this a valid thing to do?

    % clean slate
    clc; clear all; close all;

    % tspan - reversed!
    R = 0.25:-0.001:0;

    % initial values
    hinits=[0.0000001;0]; 

    % solve
    [R,v] = ode45(@equ7,R,hinits);

    % strip imaginary values (can't plot 'em)
    v(find(real(v)~=v)) = NaN;

    % plot first column 
    plot(R,v(:,1)); 

    function vprime = equ7(R,v);
    % Solve second order non-linear differential equation 7:
    %   v''(R) + 2/R*v'(R) = K_plus/(R^2)*( v^(-1/2) - lamda_plus*(1-v)^(-1/2)
    %
    % Matlab ode45 only likes first derivatives, so let:  
    %    v_1(R) = v(R)
    %    v_2(R) = v'(R)
    %
    % And create a system of first order diff eqs:
    %    v_1'(R) = v_2(R)
    %    v_2'(R) = -2/R*v_2(R) + K_plus/(R^2)*( v_1(R)^(-1/2) - lamda_plus*(1-v_1(R))^(-1/2)
    %
    % Constant Parameters:
    K_plus = 0.0859;
    lambda_plus = 3.7;

    % Build result in pieces for easier debugging of problematic terms
    int1 = 1 - v(1);
    int2 = int1^(-1/2);
    int3 = v(1)^(-1/2);
    int4 = K_plus/(R^2);
    vprime2 = -2/R*v(2);
    vprime2 = vprime2 + int4*( int3 - lambda_plus*(int2) );
    vprime=[v(2); vprime2 ];

[Jan]

A clear all removes all breakpoints of the debugger, deltes all loaded functions from the memory, such that an expensive re-import from tzhe hard disk is required. Clearing all local variables at the top aof a function is meaningless. So this command has a bunch of drawbacks and not a single benefit.

Avoid the expensive power operation int1^(-1/2) and prefer 1/sqrt(int1).

[GlobalPostcard]

Thanks for the feedback! I'll make those changes.

I'm still interested answers to my main "ode45 reversed tspan" question.

Answer 1

[Torsten]

It's a valid method to integrate backwards, and there is no loss in precision.

Usually such problems as the one above are boundary value problems with boundary conditions v'(0)=0 and v(0.25) = something. These can be solved by bvp4c.

Another way to treat your problem from above is to start integration not exactly at R=0, but at R=eps for a small number eps>0.

Best wishes

Torsten.

Answer 2

[kangkan choudhury]

Hello.

Could you post our original differential equations? Thats is without time being reversed.

The thing I wanted to understand is whether one needs to change the signs of the derivatives when one

reverses time?

Thanks and Regards

[Walter Roberson]

>> r = rand; [t,x] = ode45(@(t,y) t+y, 0:10, r)

t =

     0
     1
     2
     3
     4
     5
     6
     7
     8
     9
    10


x =

         0.814723686393179
          2.93302517351755
           10.410204360851
          32.4534055209134
          94.0925073410163
          263.366465788663
          725.227842835514
          1982.43935277454
          5401.67764062545
          14698.0253864335
           39968.110801003

>> [tr,xr] = ode45(@(t,y) t+y, fliplr(0:10), x(end))

tr =

    10
     9
     8
     7
     6
     5
     4
     3
     2
     1
     0


xr =

           39968.110801003
          14701.7124671153
          5403.14079902156
           1983.0691324623
          725.882937707868
          263.743457357354
          94.2241969828854
          32.5099185337234
           10.440162537483
          2.94490360299423
         0.819663327804449

>> r

r =

         0.814723686393179

xr(end) is close to the original boundary condition, r. Notice that we did not change anything about the equation: we just reversed time and put in the appropriate boundary condition.

[John D'Errico]

Please don't add an answer just to ask a question. But, no. There is no need to reverse the sign on the derivatives. The negative step created by a reversed tspan is all that is needed.

[kangkan choudhury]

Thanks a lot for the clarification.