Solving a system of second order ODE backwards

조회 수: 89 (최근 30일)
Padi
Padi 2023년 3월 2일
댓글: John D'Errico 2023년 3월 4일
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.)

채택된 답변

John D'Errico
John D'Errico 2023년 3월 2일
편집: John D'Errico 2023년 3월 2일
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
Padi
Padi 2023년 3월 2일
How can you impose what time step you want? I want it to be precisely h.
Thanks!
John D'Errico
John D'Errico 2023년 3월 4일
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

댓글을 달려면 로그인하십시오.

추가 답변 (3개)

Torsten
Torsten 2023년 3월 2일
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
Torsten
Torsten 2023년 3월 2일
편집: Torsten 2023년 3월 2일
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
Padi 2023년 3월 2일
I see! thanks!

댓글을 달려면 로그인하십시오.


William Rose
William Rose 2023년 3월 2일
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.
  댓글 수: 2
Padi
Padi 2023년 3월 2일
Thank you very much William!
William Rose
William Rose 2023년 3월 2일
@Padi, you're welcome.

댓글을 달려면 로그인하십시오.


Padi
Padi 2023년 3월 2일
Thank you all! This is very helpful!

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by