Solving a system of second order ODE backwards

82 visualizzazioni (ultimi 30 giorni)
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.)

Risposta accettata

John D'Errico
John D'Errico il 2 Mar 2023
Modificato: John D'Errico il 2 Mar 2023
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.
  2 Commenti
Padi
Padi il 2 Mar 2023
How can you impose what time step you want? I want it to be precisely h.
Thanks!
John D'Errico
John D'Errico il 4 Mar 2023
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

Accedi per commentare.

Più risposte (3)

Torsten
Torsten il 2 Mar 2023
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
  4 Commenti
Torsten
Torsten il 2 Mar 2023
Modificato: Torsten il 2 Mar 2023
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.

Accedi per commentare.


William Rose
William Rose il 2 Mar 2023
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
Padi il 2 Mar 2023
Thank you all! This is very helpful!

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by