# How to get the solution to a GIVEN input VECTOR using ode45 solver

2 views (last 30 days)
Ivan on 31 Jul 2020
Commented: J. Alex Lee on 3 Aug 2020
Hi everyone!
I already posted my problem here, but unfortunately we did not come to a proper solution, so I've decided to reformulate my problem in an easier way, hoping it to be clear enough.
In the code below, you'll find the variable u_opt which represents the imput u of my system of equations.
u in function calc_only_f is a single scalar number, not a vector.
u_opt is a 50x1 double vector I obtained from a previous calculation as a solution of an optimal control problem (THIS INFO IS NOT RELEVANT). What I want to acheive is using u_opt as the imput of the ode45 solver and get for EACH component of this vector (50x1 i.e. for 50 different values) a 4x1 double solution as output.
After N=50 interation I would like to have a 50x4 double solution (aka trajectory) for those 50x1 double inputs of vector u_opt
As You can see below (represented by an arrow <-- in second to last line) I obtain for EACH component (i.e. value) of the input vector u_opt a 50x4 double solution. So for N=50 interations, I get 50 of these 50x4 double trajectoies, which are given as 1x50 double cell array.
Is there possibility to acheive exactly what I'm looking for, so getting a 50x4 double solution for those 50x1 double inputs? I spent a lot of hours in finding a solution, but unfortunately I'm still stuck on this problem.
........... calcf .........
% System of equations x' = f(t,x,u)
function f = calc_only_f(t,x,u,param)
% Parameter of the systems given as a struct
M = param.M;
g = param.g;
L = param.L;
m = param.m;
d = param.d;
Sx = sin(x(3)); % sin(x)
Cx = cos(x(3)); % cos(x)
N = M + m*(1-Cx^2); % denominator
K = N*L; % L*(M + m*(1-Cx^2))
f = [x(2); (1/N)*(-m*g*Cx*Sx + m*L*x(4)^2*Sx - d*x(2)+ u); x(4); (1/(N*L))*((m+M)*g*Sx -Cx*(m*L*x(4)^2*Sx - d*x(2)) - Cx*u)];
end
% +++++++++++++++++++++++++++++++++++++++++++++++++++
% -------------------- MAIN -------------------------
% +++++++++++++++++++++++++++++++++++++++++++++++++++
N=50; sizeu=1; tf= 5;
u_opt = rand(N*sizeu,1,'double'); % only for the implementation ... dimension of u_opt should be a 50x1 double
x0 = [0;0;pi;0]; % start point 4x1 double ... x0_1 position x0_2 velocity x0_3 angle x0_4 angular velocity
param.m = 1.5;
param.M = 27;
param.L = 2;
param.g = -9.80665;
param.d = 1.2;
param.N = 50;
param.timePart = linspace(t0,tf,N);
tspan = param.timePart;
for k = 1:numel(u_opt)
[t_ode{k},x_ode{k}] = ode45(@(t,x)calc_only_f(t,x,u_opt(k),param),tspan,x0); % <-- I get a 1x50 double cell array
end
I appreciate any help!
Cheers !
Ivan on 2 Aug 2020
"It's going to give you solutions at all the time steps you requested in tSpan, which in this case is 50 long." --> Yes and as I wrote previoulsly, I got a 1x50 cell array: each cell contains a 50x4 solution for each element of t_opt. This result does not help me because so I do not have anything to compare with my reference.
Yes, you got it right!
For each single value of u drawn from u_opt I want ode45 to produce a 1x4 solution (make sure not to confuse rows and columns ;) ). So for 50 calls, as You say, having a 50x4 double solution.
That's why I'm looking for a 50x4 solution for a given 50x1 input.
Considering only one time instant (for example at a time t_c) was only a hint for the community. This would help me at least to get a trajectory at that time instant t_c, let consider t_c to be at the half of tspan, for instance.
Maybe You or someothers know how to code it.

J. Alex Lee on 2 Aug 2020
If you're just looking for code to extract some element out of your ode45 solution:
tAll = nan(1,length(u_opt));
xAll = nan(4,length(u_opt));
idx = 50; % or whatever index of tspan you want
for k = 1:numel(u_opt)
[t,x] = ode45(@(t,x)calc_only_f(t,x,u_opt(k),param),tspan,x0);
tAll(k) = t(idx);
xAll(:,k) = x(:,idx);
end
and just change idx to whatever you want. If you want t_c to be "half of tspan", and you have defined tspan with linspace, then let idx=25, for example.
J. Alex Lee on 3 Aug 2020
Ok, glad it gets you closer to what you want.
But for would-be consumers of this "solution", be warned that in my understanding this was essentially a long protracted discussion about how to extract elements from a matrix, so that is all this answer purports to provide.

R2020a

### Community Treasure Hunt

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

Start Hunting!