Graphs using ode45 do not look like they should. Why is time vector not a linear set?

조회 수: 10 (최근 30일)
I am trying to solve a system of differential equations of variables x, v, temperature with respect to time. However, when I use the ode45 function it seems like the time values that are not a linear set. When solving the system with ode45 and graphing velocity vs time, this is the graph that is output.
However, when I use the linspace function to create a linear set of data points to use as my time values, I get this graph.
The graph is sort of messed up for whatever reason, but you can see that the velocity seems to move more symmetrically. In the context of the problem that I am trying to solve, velocity should look closer to the second graph (without the weird tails on the end). I am not sure why the ode45 function is creating such a strange velocity graph. Also, when I graph x vs time, it looks like x moves symmetrically (as it should). I don't know why the velocity graph is not simply the derivative of x, as I define it within my program. Could someone help me with this problem?
My code is below
p_i = 4 ;
y0 = [1 0 1] ; %initial values
tspan = [0 40] ;
[t, y] = ode45(@rate,tspan,y0) ;
x = y(:,1) ;
v = y(:,2) ;
temp = y(:,3) ;
p = p_i * (temp ./ x) ;
%t = linspace(0,100,117) ;
figure(1)
plot(t,p)
title('pressure vs time')
figure(2)
plot(t,temp)
title('temperature vs time')
figure(3)
plot(t,v)
title('velocity vs time')
function dydt = rate(~,y)
x = y(1);
v = y(2);
temp = y(3);
p_i = 4 ; %initial p
a = 2/5 ; %constant
p = p_i * (temp ./ x) ; %pressure in termps of p_i, temp, x
dxdt = v ; %diff eq 1
dvdt = (p - 1) / (p_i - 1) ; %diff eq 2
dtempdt = -v .* (p / p_i) * a ; %diff eq 3
dydt = [dxdt ; dvdt ; dtempdt] ;
end

채택된 답변

Paul
Paul 2022년 6월 12일
편집: Paul 2022년 6월 13일
Hi Josh,
Maybe the plots are deceptive. Running the code exactly as given, except for the plotting, shows that v is the derivative of x wrt t
p_i = 4 ;
y0 = [1 0 1] ; %initial values
tspan = [0 40] ;
[t, y] = ode45(@rate,tspan,y0) ;
x = y(:,1) ;
v = y(:,2) ;
temp = y(:,3) ;
figure(1)
subplot(211);
plot(t,x)
subplot(212)
plot(t,v,t,gradient(x,t),'o')
legend('v','dx/dt')
title('velocity vs time')
function dydt = rate(~,y)
x = y(1);
v = y(2);
temp = y(3);
p_i = 4 ; %initial p
a = 2/5 ; %constant
p = p_i * (temp ./ x) ; %pressure in termps of p_i, temp, x
dxdt = v ; %diff eq 1
dvdt = (p - 1) / (p_i - 1) ; %diff eq 2
dtempdt = -v .* (p / p_i) * a ; %diff eq 3
dydt = [dxdt ; dvdt ; dtempdt] ;
end
  댓글 수: 4
Torsten
Torsten 2022년 6월 12일
Other than that, it's not clear to me what the expected result should be if not what's shown above.
As far as I understand, a symmetrical velocity field (like the position field) by a suitable change of the differential equations and/or the initial conditions.

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

추가 답변 (1개)

Walter Roberson
Walter Roberson 2022년 6월 12일
tspan = [0 40] ;
[t, y] = ode45(@rate,tspan,y0) ;
You asked ode45 to compute for some points between t = 0 and t = 40.
Hypothesize that your expectation was correct, that the times returned by ode45 are linear.
Now, construct a sync function with the pole at t = 40
x = linspace(0, 40, 250);
y = sinc(x-40) .* exp(x-40);
plot(x, y)
xlim([0 35])
figure
plot(x, y)
so pretty much flat until t = 32 or so, and then suddenly has a lot of variation.
Under your hypothesis that ode45 should return linear time... then which linear times should it return? When it processes up to 30-ish it can pretty much use points a fair distance apart. It has no idea that it needs to deal with steeper values (which requires points closer together) until it gets up past 35.
So... is your expectation that ode45() should examine the final output that it computes, and then go backwards and "fill in" samples according to the closest sampling that it turned out to need anywhere on the function, for the sake of ensuring that the output times are linear?
Or is your expectation that it will use a fixed timestep? And if so, then what fixed time-step, and what do you want ode45 to do about the locations where the values are changing rapidly where a smaller timestep might be called for?
ode45() is an adaptive ODE solver. Using variable time-steps is inherent to what it does and how it does it.
You have three options at this point:
  • switch to a fixed-timestep solver. There are some in the File Exchange, and there is a Question posted by Mathworks Support that has code for fixed-timestep solvers; OR
  • guess how far apart the important features are for your ode, and pass ode45() a vector for tspan with at least three elements, which lists the exact times that you want outputs for; OR
  • stop expecting the outputs to be at fixed timesteps, and start expecting that ode45 takes shorter steps in places where the slope is changing more quickly, and gradually starts taking longer steps again as it decides it has become safe to do so. This is inherently not symmetric in time.
  댓글 수: 2
Josh Ciesar
Josh Ciesar 2022년 6월 12일
It seems like the issue I am having is related more to the way that matlab is solving the ode system, not the time step. When I graph 'x' , the graph looks as expected. When I graph 'v', which I have defined as the velocity of x, you can see that this does not look true based on the graphs.
based on information that I have, I know that the velocity graph is supposed to oscillate symetrically. Do you know why this might be happening?
Torsten
Torsten 2022년 6월 12일
The points where velocity is zero are the points where position change its direction from up to down and vice versa. This is consistent. So I don't have a doubt that x' = v as requested.

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

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

제품


릴리스

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by