Solving a system of second order ODE backwards

Dear all,

I am trying to solve a simple second order ODE but I was hoping to solve it backwards. With that I mean, that I have info at t=1 and I want to get the value of the solution for t in (0,1]. I believe that myODE will solve it forward...meaning you need an initial condition for t=0. Could you help me out here?

function dv=myODE(h,y)
% v = v_tt + 1/t v_t+u
%where u is known and v_t=0 at t=1.
t = (0+h :h: 1-h);  
u=1+cos(2*pi*t); 
dy(1) = y(2);
dy(2)=y(1)-h*y(2)-u;
end

Also, when calling this function, I am not sure why it is giving me an error

[t,v]=ode45(@myode((t,v),h:T,[1;0]) (Here I was trying to put some intiial conditions because Im not sure how to do it backwards, but it still gives me errors forward in time.)

Answer 1

[John D'Errico]

Easy, peasy. For example, solve the ODE

y' = sin(t)

subject to the initial condition, that at t==1 we have y(1)=1/2. Now solve it moving from 1 to 0.

odefun = @(t,y) sin(t);
y0 = 1/2;
tspan = [1 0]; % it will solve with a negative time step
[tout,yout] = ode45(odefun,tspan,y0);
plot(tout,yout,'r-',1,y0,'bo')

So it started at t==1, and then used a negative time step, moving down to t==0.

Again, the initial condition is indeed y(1)==1/2. You can see the solution passes through that point.

[Padi]

How can you impose what time step you want? I want it to be precisely h.

Thanks!

[John D'Errico]

ODE45 is an adaptive solver, so your time steps are not truly the steps used by ODE45. It may need to go smaller steps, or it may interpolate as needed. But if you specify specific time steps, ODE45 will give them to you.

odefun = @(t,y) sin(t);
y0 = 1/2;
tspan = 1:-0.1:0;
[tout,yout] = ode45(odefun,tspan,y0);
[tout,yout]

ans = 11×2

1.0000 0.5000 0.9000 0.4187 0.8000 0.3436 0.7000 0.2755 0.6000 0.2150 0.5000 0.1627 0.4000 0.1192 0.3000 0.0850 0.2000 0.0602 0.1000 0.0453

Answer 2

[Torsten]

An example:

fun = @(t,y) y;
[t1,y1] = ode45(fun,0:0.01:1,1);
[t2,y2] = ode45(fun,1:-0.01:0,exp(1));
hold on
plot(t1,y1)
plot(t2,y2)
hold off
grid on

[Padi]

This is very helpful. I see how to do the precise time-step from here. Thank you very much!!

[Steven Lord]

Note that specifying a tspan vector with more than 2 elements does not control which steps the ODE solver uses internally. It specified for which times the ODE solver returns the results.

[Torsten]

@Padi

The vector "tspan" you specify (here: 0:0.01:1 or 1:-0.01:0) is only used for outputting the solution. The solver will generate its own time steps internally depending on the difficulty of the problem in certain regions of the interval of integration. So you cannot prescribe a stepsize h for the solver.

[Padi]

I see! thanks!

Answer 3

[William Rose]

@Padi,

First, try to get your script to work in the normal forward direction.

The funciton needs to be adjusted. You have

function dv=myODE(h,y)
% v = v_tt + 1/t v_t+u
%where u is known and v_t=0 at t=1.
t = (0+h :h: 1-h);  
u=1+cos(2*pi*t); 
dy(1) = y(2);
dy(2)=y(1)-h*y(2)-u;
end

That will not work for several reasons.

1. The output variable name, dv, does not match dy(1) and dy(2) used in the function body.
2. dy must be a column vector, not a row vector.
3. The internal definition of t as a vector is not needed, and causes trouble.
4. You need to pass t and y and, possibly, h, where h is a constant of the model. Alternatively, you could define h inside the function, which is what I do below. I assume u is an external forcing function.

A better version of the function would be

function dy=myODE(t,y)
% v = v_tt + 1/t v_t+u
% where u is known and v_t=0 at t=1.
dy=[0;0];           % assure dy is a column vector
h=1;                % adjust as desired
u=1+cos(2*pi*t);    % external forcing
dy(1) = y(2);
dy(2) = y(1)-h*y(2)-u;
end

Your call to ode45() includes "myode" but its name is "myODE". The case matters.

tspan=[0,1];        % time span, going forward
y0=[1,0];           % initial conditions
[t,v]=ode45(@myODE,tspan,y0);

Try that. Once it is working in the forward direction, then we can think about running it backward.

One option for going backwards is to define t2=1-t. Then reverse the signs of the derivatives in myODE (since that is what will happen when you do the substitution of variables), then integrate t2 forward, from t2=0 to 1. Specify the initial condition (at time t2=0) as whatever the state is at t=1. When you get the answer [t2,y] from ode45(), compute t=1-t2, and plot t versus y.

[Padi]

Thank you very much William!

[William Rose]

@Padi, you're welcome.

Answer 4

[Padi]

Thank you all! This is very helpful!