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

3 views (last 30 days)
Josh Ciesar
Josh Ciesar on 12 Jun 2022
Edited: Paul on 13 Jun 2022
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

Accepted Answer

Paul
Paul on 12 Jun 2022
Edited: Paul on 13 Jun 2022
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 Comments

Sign in to comment.

More Answers (1)

Walter Roberson
Walter Roberson on 12 Jun 2022
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 Comments
Torsten
Torsten on 12 Jun 2022
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.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by