Errors using ode45, state representation with time variant parameter matrices

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);

Answer 1

[Paul]

%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

  Name      Size            Bytes  Class     Attributes

  A         8x8                 8  sym                 
  p         8x1                64  double              

A has only one symbolic variable

symvar(A)

ans = 

t

Convert it to an anonymous function that takes in double t and returns A(t) as a double using matlabFunction

Afunc = matlabFunction(A)

Afunc = function_handle with value:
    @(t)reshape([0.0,0.0,0.0,0.0,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^4.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^4.*(-6.496875e+4)-(t.*(9.0./2.0e+2)+9.0./4.0e+2).^5.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^5.*2.8153125e+5-t.^4.*(t.*(9.0./2.0e+2)-1.0).^4.*2.6641248046875e-1-t.^5.*(t.*(9.0./2.0e+2)-1.0).^5.*5.195043369140625e-2-8.12109375e-2,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^4.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^4.*2.8153125e+5+(t.*(9.0./2.0e+2)+9.0./4.0e+2).^5.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^5.*1.2387375e+6+t.^4.*(t.*(9.0./2.0e+2)-1.0).^4.*1.15445408203125+t.^5.*(t.*(9.0./2.0e+2)-1.0).^5.*2.285819082421875e-1+2.8153125e-1,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^2.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^2.*2.0e+2+t.^2.*(t.*(9.0./2.0e+2)-1.0).^2.*(8.1e+1./2.0e+2),(t.*(9.0./2.0e+2)+9.0./4.0e+2).^2.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^2.*7.03125e+3-t.^2.*(t.*(9.0./2.0e+2)-1.0).^2.*1.423828125e+1,0.0,0.0,0.0,0.0,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^5.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^5.*6.496875e+4+(t.*(9.0./2.0e+2)+9.0./4.0e+2).^6.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^6.*2.8153125e+5+t.^5.*(t.*(9.0./2.0e+2)-1.0).^5.*1.198856162109375e-2+t.^6.*(t.*(9.0./2.0e+2)-1.0).^6.*2.337769516113281e-3+1.0828125e-1,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^5.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^5.*(-2.8153125e+5)-(t.*(9.0./2.0e+2)+9.0./4.0e+2).^6.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^6.*1.2387375e+6-t.^5.*(t.*(9.0./2.0e+2)-1.0).^5.*5.195043369140625e-2-t.^6.*(t.*(9.0./2.0e+2)-1.0).^6.*1.028618587089844e-2-4.926796875e-1,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^3.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^3.*-2.0e+2-t.^3.*(t.*(9.0./2.0e+2)-1.0).^3.*1.8225e-2,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^3.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^3.*(-7.03125e+3)+t.^3.*(t.*(9.0./2.0e+2)-1.0).^3.*6.4072265625e-1,0.0,0.0,0.0,0.0,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^2.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^2.*6.496875e+4+(t.*(9.0./2.0e+2)+9.0./4.0e+2).^3.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^3.*2.8153125e+5+t.^2.*(t.*(9.0./2.0e+2)-1.0).^2.*1.3156171875e+2+t.^3.*(t.*(9.0./2.0e+2)-1.0).^3.*2.565453515625e+1,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^2.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^2.*(-2.8153125e+5)-(t.*(9.0./2.0e+2)+9.0./4.0e+2).^3.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^3.*1.2387375e+6-t.^2.*(t.*(9.0./2.0e+2)-1.0).^2.*5.7010078125e+2-t.^3.*(t.*(9.0./2.0e+2)-1.0).^3.*1.128799546875e+2,-4.0e+2,0.0,0.0,0.0,0.0,0.0,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^2.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^2.*2.92359375e+5+(t.*(9.0./2.0e+2)+9.0./4.0e+2).^3.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^3.*1.266890625e+6-t.^2.*(t.*(9.0./2.0e+2)-1.0).^2.*5.92027734375e+2-t.^3.*(t.*(9.0./2.0e+2)-1.0).^3.*1.15445408203125e+2,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^2.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^2.*(-1.266890625e+6)-(t.*(9.0./2.0e+2)+9.0./4.0e+2).^3.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^3.*5.57431875e+6+t.^2.*(t.*(9.0./2.0e+2)-1.0).^2.*2.565453515625e+3+t.^3.*(t.*(9.0./2.0e+2)-1.0).^3.*5.0795979609375e+2,0.0,-6.328125e+4,1.0,0.0,0.0,0.0,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^4.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^4.*(-4.67775e+2)-(t.*(9.0./2.0e+2)+9.0./4.0e+2).^5.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^5.*2.027025e+3-t.^4.*(t.*(9.0./2.0e+2)-1.0).^4.*1.918169859375e-3-t.^5.*(t.*(9.0./2.0e+2)-1.0).^5.*3.74043122578125e-4,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^4.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^4.*2.027025e+3+(t.*(9.0./2.0e+2)+9.0./4.0e+2).^5.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^5.*8.91891e+3+t.^4.*(t.*(9.0./2.0e+2)-1.0).^4.*8.312069390625e-3+t.^5.*(t.*(9.0./2.0e+2)-1.0).^5.*1.64578973934375e-3,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^2.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^2.*(3.6e+1./2.5e+1)+t.^2.*(t.*(9.0./2.0e+2)-1.0).^2.*2.916e-3,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^2.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^2.*(4.05e+2./8.0)-t.^2.*(t.*(9.0./2.0e+2)-1.0).^2.*1.02515625e-1,0.0,1.0,0.0,0.0,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^5.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^5.*4.67775e+2+(t.*(9.0./2.0e+2)+9.0./4.0e+2).^6.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^6.*2.027025e+3+t.^5.*(t.*(9.0./2.0e+2)-1.0).^5.*8.6317643671875e-5+t.^6.*(t.*(9.0./2.0e+2)-1.0).^6.*1.683194051601563e-5,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^5.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^5.*(-2.027025e+3)-(t.*(9.0./2.0e+2)+9.0./4.0e+2).^6.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^6.*8.91891e+3-t.^5.*(t.*(9.0./2.0e+2)-1.0).^5.*3.74043122578125e-4-t.^6.*(t.*(9.0./2.0e+2)-1.0).^6.*7.406053827046875e-5,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^3.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^3.*(-3.6e+1./2.5e+1)-t.^3.*(t.*(9.0./2.0e+2)-1.0).^3.*1.3122e-4,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^3.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^3.*(-4.05e+2./8.0)+t.^3.*(t.*(9.0./2.0e+2)-1.0).^3.*4.613203125e-3,0.0,0.0,1.0,0.0,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^2.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^2.*4.67775e+2+(t.*(9.0./2.0e+2)+9.0./4.0e+2).^3.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^3.*2.027025e+3+t.^2.*(t.*(9.0./2.0e+2)-1.0).^2.*9.47244375e-1+t.^3.*(t.*(9.0./2.0e+2)-1.0).^3.*1.84712653125e-1,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^2.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^2.*(-2.027025e+3)-(t.*(9.0./2.0e+2)+9.0./4.0e+2).^3.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^3.*8.91891e+3-t.^2.*(t.*(9.0./2.0e+2)-1.0).^2.*4.104725625-t.^3.*(t.*(9.0./2.0e+2)-1.0).^3.*8.1273567375e-1,-7.2e+1./2.5e+1,0.0,0.0,0.0,0.0,1.0,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^2.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^2.*2.1049875e+3+(t.*(9.0./2.0e+2)+9.0./4.0e+2).^3.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^3.*9.1216125e+3-t.^2.*(t.*(9.0./2.0e+2)-1.0).^2.*4.2625996875-t.^3.*(t.*(9.0./2.0e+2)-1.0).^3.*8.312069390625e-1,(t.*(9.0./2.0e+2)+9.0./4.0e+2).^2.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^2.*(-9.1216125e+3)-(t.*(9.0./2.0e+2)+9.0./4.0e+2).^3.*(t.*(9.0./2.0e+2)-3.91e+2./4.0e+2).^3.*4.0135095e+4+t.^2.*(t.*(9.0./2.0e+2)-1.0).^2.*1.84712653125e+1+t.^3.*(t.*(9.0./2.0e+2)-1.0).^3.*3.657310531875,0.0,-4.55625e+2],[8,8])

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)

[Michael]

Thank you so much for your help!

[Paul]

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