How do I get I pass out the input parameter to ode45 integration

2 views (last 30 days)
I want to integrate velocity and obtain posiiton. Is there a way to do this with ode45 and obtain the set of velocities as additonal output to the ode45. Here is a code I have written but I am not sure it works correctly.
close all; clear all;clf
x0=[0.5;0.5];
tf=5; N=10000
tspan = linspace(0, tf,N);
dt =tf/N;
[t,x]= ode45(@(t,x)position_calc(x), tspan, x0);
figure
plot(t,x(:,1),'b--','linewidth',2)
hold on
% figure
plot(t,x(:,2),'g--');
grid
legend(["x","vel","pos"])
function [dstate]= position_calc(state)
x=state(1);
dx= state(2);
vel=x^2;
dstate= [x;vel];
end
I do not want to run the function inside ode again to obtain the velocity
  3 Comments
Hafeez Jimoh
Hafeez Jimoh on 31 Jan 2023
vel= x^2. But I also need the values of x itself. I want to integrate velocity to obtain position. So x(:,2) is intended to be velocity and x(:,1) is intended to be position. Kindly let me know if I am doing something wrong. Thanks
Torsten
Torsten on 31 Jan 2023
Edited: Torsten on 31 Jan 2023
If x(:,1) is intended to be position and x(:,2) is intended to be velocity, your odes read
y''(t) = a
where "a" is acceleration.
Written as a system
y1' = y2
y2' = a
where y1 is position and y2 is velocity.
So I don't understand your function "position_calc".

Sign in to comment.

Accepted Answer

William Rose
William Rose on 31 Jan 2023
I think @Torsten, is right, as usual. I think you are simulating a one-dimensional system in which the velocity equals the square of the position:
v = x^2.
Therefore I recommend that x represent position only. In other words, x should be a vector, rather than an array.
After you do the integration, you can compute the velocity.
Here is a code fragment that integrates the system for 1.5 seconds. I find that by t=2, it blows up.
x0=0.5; %initial position
tf=1.5; %final time
[t,x]= ode45(@(t,x) x*x, [0 tf], x0);
v=x.*x; %velocity
plot(t,x,'b',t,v,'-r') %plot results
grid on; legend(["Position","Velocity"]); xlabel('Time')
Good luck.
  3 Comments
Sam Chak
Sam Chak on 31 Jan 2023
Edited: Sam Chak on 31 Jan 2023
Hafeez, I suggest you show the Differential equation that you want to solve. Only math, no MATLAB code. This will allow us to clearly look into the math problem and then check if it can be solved using ode45 solver.
What you described sounds like a kinematic differential equation. But what you want to obtain the velocity state directly from the ode45 solution indicates a dynamical differential equation.
Torsten
Torsten on 31 Jan 2023
Is there anyway to avoid doing this :
v=x.*x; %velocity
and getting the velocity directly from the ode45 line?
Not really. The usual way is to call the function in which you supply the derivatives with the result vectors from ode45 after integration has finished:
x0 = 0.5; %initial position
tf = 1.5; %final time
[t,x] = ode45(@fun, [0 tf], x0);
v = fun(t,x); %velocity
plot(t,x,'b',t,v,'-r') %plot results
grid on; legend(["Position","Velocity"]); xlabel('Time')
function dx = fun(t,x)
dx = x.^2;
end

Sign in to comment.

More Answers (1)

Sulaymon Eshkabilov
Sulaymon Eshkabilov on 31 Jan 2023
There are a couple of overlooked points as Torsten pinpointed out. Here is one simple explae of dx = 2x example solution very similar to your exercise:
% Let's say: dx = 2*x; which is the velocity formulation: f(t,x) = 2*x;
% then solution x(t) can be found using dsolve() that gives an analytical
% solution or ode45 gives numerical solution
close all; clearvars
x0=0.5; % Note: dx has one Initial Condition, i.e., x(0) = 0.5, for example
tf=5; N=10000;
tspan = linspace(0, tf,N);
[t,x]= ode45(@(t,x)position_calc(x), tspan, x0);
figure
yyaxis left
plot(t,x,'bo-','linewidth',2, 'DisplayName','Displacement')
ylabel('$x(t), \ [m]$', 'Interpreter','latex')
% Now, the velocity values can be obtained from the integrated displacement values,
v =diff(x);
yyaxis right
plot(t(1:end-1),v,'r-.', 'LineWidth',1.75, 'DisplayName','Velocity');
ylabel('$v(t), \ [m/s]$', 'Interpreter','latex')
xlabel('$time, \ [s]$', 'Interpreter','latex')
grid on
legend('Show', 'Location', 'Best')
figure
D = array2table(t(1:end-1));
D = renamevars(D, 'Var1', 'time');
D.x = x(1:end-1);
D.v = v;
H=stackedplot(D, 'XVariable', 'time');
H.LineWidth=2;
H.LineProperties(1).Color = 'b';
H.LineProperties(1).Marker = 'o';
H.LineProperties(2).Color = 'r';
H.LineProperties(2).LineStyle = '-.';
grid on
function dx= position_calc(state)
dx= 2*state;
end

Community Treasure Hunt

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

Start Hunting!