Errors using ode45, state representation with time variant parameter matrices
Show older comments
I am trying to model a car moving across a beam using the assumed modes method and I am getting an error with ode45. For the system, the admissible functions are time varying and go into the parameter matrices that form the state representation that is being solved for. To include variable t in the admissible functions, as opposed to trying to construct all of the matrices in the ode45 command line, I used syms to add in the variable, not sure if that is the correct way to do so or if that would be causing the errors. When I run the script, I get the following errors:
Error using superiorfloat
Inputs must be floats, namely single or double.
Error in odearguments (line 114)
dataType = superiorfloat(t0,y0,f0);
Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
Error in Proejct_2B_Script (line 53)
[t,z] = ode45(@(t,z) A*z + p, tspan, iniCon);
Here is the script:
%System Parameters
a = 4.5; M = 2.5e3; I = 3.2e2;
k = 5e5; c = 3.6e3; v0 = 18; g = 9.81;
rho = 8e2; EI = 7e8; L = 400;
%Parameter Matrices
%Beam
Mm = rho*L*[1/630 1/2772; 1/2772 1/12012];
Km = EI/(L^3)*[4/5 6/35; 6/35 2/35];
%Rigid Body
Mr = [M 0; 0 I];
Dr = c*[1 -a; -a a^2]+c*[1 a; a a^2];
Kr = k*[1 -a; -a a^2]+k*[1 a; a a^2];
%Admissible function inputs
syms t
Phix1 = [((v0*t/L)^2)*(1-(v0*t/L))^2 ((v0*t/L)^3)*(1-(v0*t/L))^3];
Phix1T = [((v0*t/L)^2)*(1-(v0*t/L))^2; ((v0*t/L)^3)*(1-(v0*t/L))^3];
Phix2 = [(((v0*t+2*a)/L)^2)*(1-((v0*t+2*a)/L))^2 (((v0*t+2*a)/L)^3)*(1-((v0*t+2*a)/L))^3];
Phix2T = [(((v0*t+2*a)/L)^2)*(1-((v0*t+2*a)/L))^2; (((v0*t+2*a)/L)^3)*(1-((v0*t+2*a)/L))^3];
%Combined
Dc = c*Phix1T*Phix1+c*Phix2T*Phix2;
Dbr = c*Phix1T*[1 -a]+c*Phix2T*[1 a];
Drb = transpose(Dbr);
Kc = k*Phix1T*Phix1+k*Phix2T*Phix2;
Kbr = k*Phix1T*[1 -a]+k*Phix2T*[1 a];
Krb = transpose(Kbr);
%Assembly Matrices
Ma = [Mm zeros(2,2); zeros(2,2) Mr];
Ca = [Dc -Dbr; -Drb Dr];
Ka = [Km+Kc -Kbr; -Krb Kr];
Fa = [0; 0; -M*g; 0];
%State Representation Matrices
Minv = Ma^-1;
p = [0; 0; 0; 0; Minv*Fa];
Ident = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1];
A = [zeros(4,4) Ident; -Minv*Ka -Minv*Ca];
% %ODE Solver
tspan = [0 (L-2*a)/v0];
iniCon = [0; 0; .035; 0.015/a; 0; 0; 0.15; -0.45/a];
[t,z] = ode45(@(t,z) A*z + p, tspan, iniCon);
Accepted Answer
%System Parameters
a = 4.5; M = 2.5e3; I = 3.2e2;
k = 5e5; c = 3.6e3; v0 = 18; g = 9.81;
rho = 8e2; EI = 7e8; L = 400;
%Parameter Matrices
%Beam
Mm = rho*L*[1/630 1/2772; 1/2772 1/12012];
Km = EI/(L^3)*[4/5 6/35; 6/35 2/35];
%Rigid Body
Mr = [M 0; 0 I];
Dr = c*[1 -a; -a a^2]+c*[1 a; a a^2];
Kr = k*[1 -a; -a a^2]+k*[1 a; a a^2];
%Admissible function inputs
syms t
Phix1 = [((v0*t/L)^2)*(1-(v0*t/L))^2 ((v0*t/L)^3)*(1-(v0*t/L))^3];
Phix1T = [((v0*t/L)^2)*(1-(v0*t/L))^2; ((v0*t/L)^3)*(1-(v0*t/L))^3];
Phix2 = [(((v0*t+2*a)/L)^2)*(1-((v0*t+2*a)/L))^2 (((v0*t+2*a)/L)^3)*(1-((v0*t+2*a)/L))^3];
Phix2T = [(((v0*t+2*a)/L)^2)*(1-((v0*t+2*a)/L))^2; (((v0*t+2*a)/L)^3)*(1-((v0*t+2*a)/L))^3];
%Combined
Dc = c*Phix1T*Phix1+c*Phix2T*Phix2;
Dbr = c*Phix1T*[1 -a]+c*Phix2T*[1 a];
Drb = transpose(Dbr);
Kc = k*Phix1T*Phix1+k*Phix2T*Phix2;
Kbr = k*Phix1T*[1 -a]+k*Phix2T*[1 a];
Krb = transpose(Kbr);
%Assembly Matrices
Ma = [Mm zeros(2,2); zeros(2,2) Mr];
Ca = [Dc -Dbr; -Drb Dr];
Ka = [Km+Kc -Kbr; -Krb Kr];
Fa = [0; 0; -M*g; 0];
%State Representation Matrices
Minv = Ma^-1;
p = [0; 0; 0; 0; Minv*Fa];
Ident = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1];
A = [zeros(4,4) Ident; -Minv*Ka -Minv*Ca];
At this point, A is a sym and p is double.
whos p A
A has only one symbolic variable
symvar(A)
Convert it to an anonymous function that takes in double t and returns A(t) as a double using matlabFunction
Afunc = matlabFunction(A)
Use Afunc in the odefun and now ode45 runs without error. You might also want to look into using simplify on A before creating Afund or on the input to matlabFunction.
% %ODE Solver
tspan = [0 (L-2*a)/v0];
iniCon = [0; 0; .035; 0.015/a; 0; 0; 0.15; -0.45/a];
[t,z] = ode45(@(t,z) Afunc(t)*z + p, tspan, iniCon);
plot(t,z)
2 Comments
Michael
on 6 Dec 2023
Paul
on 6 Dec 2023
You're very welcome.
An alternative approach would be to just put all the computations into a single function that computes zdot numerically as a function of t and z
[t,z] = ode45(@odefun,tspan, iniCon);
function zdot = odefun(t,z)
% put all of the computations here for A and p, same as above but don't use
% syms t
% compute zdot
zdot = A*z + p
end
More Answers (0)
Categories
Find more on Numeric Solvers in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
