Implementation issues in function with symbolic ode using ode45 Name returns a vector of length 1 but the length of initial conditions vector is 12.

I've been circling around errors all afternoon and after pouring over the mathwork forums and answers section I've finally settled on these couple of errors that I'm unsure of how to solve.

The idea is that I have some system of equations which have components that vary with time and I'm symbolically constructing the equations before returning this system back to my main function running the ode45 solver. I keep oscillating between errors like returns vector of length 1 but the length of intial conditions vector is 12 when I use the odefun = matlabFunction seen below in the code and Error using superiorfloat, Inputs must be floats, namely single or double when I just use [dydt,Y] = odeToVectorField. I'm not sure how to package my return from my function properly to keep all the functions happy :/

The main function calling

function [T,Y] = call_osc_1BodyKane2Verification011823()
tspan = [0 10];
g = 9.81; %m/s2
initVec = [0 0 0 0 0 0 100 -g 0 0 0 0];
[T,A] = ode45(@osc_1BodyKane_2Verification011823,tspan,initVec);
%[T,A] = ode45(@(T,A) osc_1BodyKane_2Verification011823(T,A),tspan,initVec);

(Plotting stuff here)

end

and then the function I'm calling

function odefun = osc_1BodyKane_2Verification011823(time,Y)
%function dydt = osc_1BodyKane_2Verification011823(time,Y)
%function [dydt,Y] = osc_1BodyKane_2Verification011823(time,Y)

syms q1(t) q2(t) q3(t) q4(t) q5(t) q6(t)

%generalized speeds 
%(q1 = xi, q2 = yi, q3 = zi, q4 = theta, q5 = phi, q6 = psi)
%dydt = zeros(12,1);

g = 9.81; %m/s^2 
rho = 1.225; %kg/m^3
mB = 10; %Generic
mPlate = 6;
rPlate = 6;
hPlate = 1;

%Calculate Forces
FLaunch = getLaunchForce(time);
FAero = getAeroForce(Y,rho)
FWeight = mB*g
FThrust = FAero(1)
Fbx = FLaunch + FAero(1)+FThrust;
Fby = FAero(2);
Fbz = FAero(3) + -FWeight;
% Fbx = 0;
% Fby = 0;
% Fbz = 0;

%Mass Moments
Ib = zeros(3,3);
IBxx = 1/12*mB*(3*rPlate*rPlate + hPlate*hPlate)
IByy = 1/12*mPlate*(3*rPlate*rPlate + hPlate*hPlate)
IBzz = 1/2*mPlate*rPlate*rPlate
Lb = 0;
Mb = 0;
Nb = 0;

u4 = diff(q5) - diff(q6)*sin(q4);
u5 = diff(q4)*cos(q5) + diff(q6)*cos(q4)*sin(q5);
u6 = diff(q6)*cos(q4)*cos(q5) - diff(q4)*sin(q5);
u1 = diff(q1) + u5*q3 - u6*q2;
u2 = diff(q2) - u4*q3 + u6*q1;
u3 = diff(q3) + u4*q2 - u5*q1;

u1dot = diff(u1);
u2dot = diff(u2);
u3dot = diff(u3);
u4dot = diff(u4);
u5dot = diff(u5);
u6dot = diff(u6);

Kanes1 = Fbx - mB*(u1dot - u2*u6 + u3*u5);
Kanes2 = Fby - mB*(u2dot + u1*u6 - u3*u4);
Kanes3 = Fbz - mB*(u3dot - u1*u5 + u2*u4);
Kanes4 = Lb + Fby*q3 - Fbz*q2 - IBxx*u4dot + IByy*u5*u6 - IBzz*u5*u6 + mB*q2*u3dot - mB*q3*u2dot - mB*q2*u1*u5 + mB*q2*u2*u4 - mB*q3*u1*u6 + mB*q3*u3*u4;
Kanes5 = Mb - 0.5000*Fbx*q3 - 0.5000*Fbx*q6 + Fbz*q1 - IByy*u5dot - 0.5000*Fbx*q3*cos(2*q5) + 0.5000*Fbx*q6*cos(2*q5) - IBxx*u4*u6 + IBzz*u4*u6 - mB*q1*u3dot + 0.5000*mB*q3*u1dot + 0.5000*mB*q6*u1dot + mB*q1*u1*u5 - mB*q1*u2*u4 - 0.5000*mB*q3*u2*u6 + 0.5000*mB*q3*u3*u5 - 0.5000*mB*q6*u2*u6 + 0.5000*mB*q6*u3*u5 + 0.5000*mB*q3*u1dot*cos(2*q5) - 0.5000*mB*q6*u1dot*cos(2*q5) - 0.5000*mB*q3*u2*u6*cos(2*q5) + 0.5000*mB*q3*u3*u5*cos(2*q5) + 0.5000*mB*q6*u2*u6*cos(2*q5) - 0.5000*mB*q6*u3*u5*cos(2*q5);
Kanes6 = Nb + Fbx*q2 - Fby*q1 - IBzz*u6dot + 0.5000*Fbx*q3*sin(2*q5) - 0.5000*Fbx*q6*sin(2*q5) + IBxx*u4*u5 - IByy*u4*u5 + mB*q1*u2dot - mB*q2*u1dot + mB*q1*u1*u6 - mB*q1*u3*u4 + mB*q2*u2*u6 - mB*q2*u3*u5 - 0.5000*mB*q3*u1dot*sin(2*q5) + 0.5000*mB*q6*u1dot*sin(2*q5) + 0.5000*mB*q3*u2*u6*sin(2*q5) - 0.5000*mB*q3*u3*u5*sin(2*q5) - 0.5000*mB*q6*u2*u6*sin(2*q5) + 0.5000*mB*q6*u3*u5*sin(2*q5);

%Y = [q4, Dq4, q1, Dq1, q6, Dq6, q3, Dq3, q5, Dq5, q2, Dq2]
[dydt,Y] = odeToVectorField(Kanes1 == 0,Kanes2 == 0, Kanes3 == 0, Kanes4 == 0, Kanes5 == 0, Kanes6 == 0)
%dydt = odeFunction(dydt,Y)
%class(dydt)
%class(Y)

odefun = matlabFunction(dydt,'vars', {t,'Y'})

end

I call the same functions, getLaunchForce and getAeroForce in other codes so I'm confindent they aren't the issue.

Thanks!

[Star Strider]

‘I call the same functions, getLaunchForce and getAeroForce in other codes so I'm confindent they aren't the issue.’

Perhaps, however we need them if we’re to run your code. We can’t run it without them.

[Megan Frisbey]

Sure

function AeroForce = getAeroForce(y,rho)
g = 9.98;
aoa = y(11)
S = 9; %Note 1 ft mac, 9 ft span
u = y(1);
v = y(2);
w = y(3);
vel = sqrt(u*u + v*v + w*w)
aero = CLCD(aoa)
CL = aero(1)
CD = aero(2)
xForce = - 0.5*CD*rho*vel*vel*S;
yForce = 0;
zForce = 0.5*CL*rho*vel*vel*S;
AeroForce = [xForce yForce zForce]
end
function LiftDrag = CLCD(aoa)
alphaL = -5;
alphaH = 10; %Linear region - From classic Clark Y for placeholder
CL_L = -0.2;
CL_H = 1.35;
profileDrag = 0.015;
oswaldE = 0.9;
AR = 9;
coeffDHEnd = sin(2*(18-45)/180*pi) +1.0;
coeffLDEnd = cos(2.0*((18-45)/180*pi));
lastLinearCD = 0;
if aoa <= alphaH && aoa > alphaL
m = (CL_H-CL_L)/(alphaH - alphaL);
b = CL_L - m*alphaL;
coeffLift = m*aoa + b;
coeffDrag = profileDrag + coeffLift^2/(pi*oswaldE*AR);
lastLinearCD = coeffDrag;
elseif aoa > alphaH && aoa <= (alphaH+8)
m = (CL_H- coeffLDEnd)/(alphaH - (8+alphaH));
coeffLift = m*aoa + CL_H - m*alphaH;
coeffDrag = lastLinearCD + (lastLinearCD-coeffDHEnd)/(alphaH - (alphaH+8))*(aoa-10);
else
coeffLift = cos(2.0*((aoa/180*pi)-45/180*pi));
coeffDrag = sin(2*((aoa-45)/180*pi))+1.0;
end
LiftDrag(1) = coeffLift;
LiftDrag(2) = coeffDrag;
end

function launchForce = getLaunchForce(t)
g = 9.98;
tLimit = 1;
if t < tLimit
launchForce = 20*g - 20*g*t
else
launchForce = 0
end
end

[Steven Lord]

You've shown a lot of code. Is there any way you could show the mathematical form of the differential equations you're trying to solve? I suspect from something you wrote in one of your comments that you don't actually have a set of ordinary differential equations but instead have a set of delay differential equations. The ODE solvers like ode45 may not be the right tool. You may need to use a DDE solver like dde23.

[Megan Frisbey]

Thank you! Yes I've wondered this also but when I hand code in the output from odeToVectorField and remove the symbolic definitions it seems to run okay. But as I'll be working with increasingly more complicated sets of equations this isn't a long term solution. I've attached the formal write up if that helps or I can go into more detail in some notes. Really what I'm constructing is equation 7, made up of equation 8 & 9 and then the relationship from the Y matrix between generalized speeds and generalized coordinates given in equation 19

I really appreciate the help!!

[Megan Frisbey]

Also just to note the dynamics equations from Kane's should produce a set of ODEs. I probably have implementation errors (I think I need to be setting 6 of my equations to equal my udot terms instead of allowing the substitution into the Kanes equations and then asking odeToVectorField to reduce them back down) which is what I've recently wondered about, but odeToVectorField yeilds a set of ODEs. If that helps explain it better.

Answer 1

[Cris LaPierre]

I would suggest moving ode45 inside your function. Here's an example borrowed from the odeToVectorField documentation page.

tspan = [0 20];
initVec = [2 0];

ySol = osc_1BodyKane_2Verification011823(tspan,initVec);

tValues = linspace(0,20,100);
yValues = deval(ySol,tValues,1);
plot(tValues,yValues)

function ySol = osc_1BodyKane_2Verification011823(tspan,initVec)
syms y(t)

eqn = diff(y,2) == (1-y^2)*diff(y)-y;

V = odeToVectorField(eqn);
M = matlabFunction(V,'vars',{'t','Y'});

ySol = ode45(M,tspan,initVec);

end

[Cris LaPierre]

The other option is to remove the '@' before the odefun input to ode45. Using the same example

tspan = [0 20];
initVec = [2 0];

[T,A] = ode45(osc_1BodyKane_2Verification011823,tspan,initVec);

plot(T,A(:,1))

function odefun = osc_1BodyKane_2Verification011823(t,Y)
syms y(t)

eqn = diff(y,2) == (1-y^2)*diff(y)-y;

V = odeToVectorField(eqn);
odefun = matlabFunction(V,'vars',{'t','Y'});

end

[Megan Frisbey]

I can move it inside and have done that previously but that doesn't allow me the option to derive symbolically a different set of ode's at each time step since the functions for getAeroForce and getLaunchForce are dependant on time and the solution from the previous time step. (unless I'm just missing something?)

When simply removing the @ I get the error Not enough input arguments. do you mean to say instead

[T,A] = ode45(osc_1BodyKane_2Verification011823(T,A),tspan,initVec);

?

Thanks

Answer 2

[Cris LaPierre]

Thank you for sharing all your functions. If your equation is changing every timestep, I feel like you can't just return the symbolic function from your odefun. I believe you need to evaluate it and return a vector of the 12 values at that time step. I think that is why you are getting an error about a 1-element vector (just a symbolic function), when a vector the same length as initVec is expected.

Modify your odefun function to evaluate odefun at the current time and with the previous values of Y using the following code:

yNew = odefun(time,Y);

Have yNew be what the function returns instead of odefun. Those outputs become the inputs for the next iteration. Time automatically increases.

This does take a long time to run (about an hour for me). You might be able to speed it up a little by adding semicolons to all your code to suppress outputs.

tspan = [0 10];
g = 9.81; %m/s2
initVec = [0 0 0 0 0 0 100 -g 0 0 0 0];
[T,A] = ode45(@osc_1BodyKane_2Verification011823,tspan,initVec);

plot(T,A(:,1))

function Ynew = osc_1BodyKane_2Verification011823(time,Y)
%function dydt = osc_1BodyKane_2Verification011823(time,Y)
%function [dydt,Y] = osc_1BodyKane_2Verification011823(time,Y)

syms q1(t) q2(t) q3(t) q4(t) q5(t) q6(t)

g = 9.81; %m/s^2
rho = 1.225; %kg/m^3
mB = 10; %Generic
mPlate = 6;
rPlate = 6;
hPlate = 1;

%Calculate Forces
FLaunch = getLaunchForce(time);
FAero = getAeroForce(Y,rho);
FWeight = mB*g;
FThrust = FAero(1);
Fbx = FLaunch + FAero(1)+FThrust;
Fby = FAero(2);
Fbz = FAero(3) + -FWeight;

%Mass Moments
Ib = zeros(3,3);
IBxx = 1/12*mB*(3*rPlate*rPlate + hPlate*hPlate);
IByy = 1/12*mPlate*(3*rPlate*rPlate + hPlate*hPlate);
IBzz = 1/2*mPlate*rPlate*rPlate;
Lb = 0;
Mb = 0;
Nb = 0;

u4 = diff(q5) - diff(q6)*sin(q4);
u5 = diff(q4)*cos(q5) + diff(q6)*cos(q4)*sin(q5);
u6 = diff(q6)*cos(q4)*cos(q5) - diff(q4)*sin(q5);
u1 = diff(q1) + u5*q3 - u6*q2;
u2 = diff(q2) - u4*q3 + u6*q1;
u3 = diff(q3) + u4*q2 - u5*q1;

u1dot = diff(u1);
u2dot = diff(u2);
u3dot = diff(u3);
u4dot = diff(u4);
u5dot = diff(u5);
u6dot = diff(u6);

Kanes1 = Fbx - mB*(u1dot - u2*u6 + u3*u5);
Kanes2 = Fby - mB*(u2dot + u1*u6 - u3*u4);
Kanes3 = Fbz - mB*(u3dot - u1*u5 + u2*u4);
Kanes4 = Lb + Fby*q3 - Fbz*q2 - IBxx*u4dot + IByy*u5*u6 - IBzz*u5*u6 + mB*q2*u3dot - mB*q3*u2dot - mB*q2*u1*u5 + mB*q2*u2*u4 - mB*q3*u1*u6 + mB*q3*u3*u4;
Kanes5 = Mb - 0.5000*Fbx*q3 - 0.5000*Fbx*q6 + Fbz*q1 - IByy*u5dot - 0.5000*Fbx*q3*cos(2*q5) + 0.5000*Fbx*q6*cos(2*q5) - IBxx*u4*u6 + IBzz*u4*u6 - mB*q1*u3dot + 0.5000*mB*q3*u1dot + 0.5000*mB*q6*u1dot + mB*q1*u1*u5 - mB*q1*u2*u4 - 0.5000*mB*q3*u2*u6 + 0.5000*mB*q3*u3*u5 - 0.5000*mB*q6*u2*u6 + 0.5000*mB*q6*u3*u5 + 0.5000*mB*q3*u1dot*cos(2*q5) - 0.5000*mB*q6*u1dot*cos(2*q5) - 0.5000*mB*q3*u2*u6*cos(2*q5) + 0.5000*mB*q3*u3*u5*cos(2*q5) + 0.5000*mB*q6*u2*u6*cos(2*q5) - 0.5000*mB*q6*u3*u5*cos(2*q5);
Kanes6 = Nb + Fbx*q2 - Fby*q1 - IBzz*u6dot + 0.5000*Fbx*q3*sin(2*q5) - 0.5000*Fbx*q6*sin(2*q5) + IBxx*u4*u5 - IByy*u4*u5 + mB*q1*u2dot - mB*q2*u1dot + mB*q1*u1*u6 - mB*q1*u3*u4 + mB*q2*u2*u6 - mB*q2*u3*u5 - 0.5000*mB*q3*u1dot*sin(2*q5) + 0.5000*mB*q6*u1dot*sin(2*q5) + 0.5000*mB*q3*u2*u6*sin(2*q5) - 0.5000*mB*q3*u3*u5*sin(2*q5) - 0.5000*mB*q6*u2*u6*sin(2*q5) + 0.5000*mB*q6*u3*u5*sin(2*q5);

%Y = [q4, Dq4, q1, Dq1, q6, Dq6, q3, Dq3, q5, Dq5, q2, Dq2]
dydt = odeToVectorField(Kanes1 == 0,Kanes2 == 0, Kanes3 == 0, Kanes4 == 0, Kanes5 == 0, Kanes6 == 0);

odefun = matlabFunction(dydt,'vars', {t,'Y'});
Ynew = odefun(time,Y); %################  Add this line to evaluate odefun at current time step

end

function AeroForce = getAeroForce(y,rho)
g = 9.98;
aoa = y(11);
S = 9; %Note 1 ft mac, 9 ft span
u = y(1);
v = y(2);
w = y(3);
vel = sqrt(u*u + v*v + w*w);
aero = CLCD(aoa);
CL = aero(1);
CD = aero(2);
xForce = - 0.5*CD*rho*vel*vel*S;
yForce = 0;
zForce = 0.5*CL*rho*vel*vel*S;
AeroForce = [xForce yForce zForce];
end

function LiftDrag = CLCD(aoa)
alphaL = -5;
alphaH = 10; %Linear region - From classic Clark Y for placeholder
CL_L = -0.2;
CL_H = 1.35;
profileDrag = 0.015;
oswaldE = 0.9;
AR = 9;
coeffDHEnd = sin(2*(18-45)/180*pi) +1.0;
coeffLDEnd = cos(2.0*((18-45)/180*pi));
lastLinearCD = 0;
if aoa <= alphaH && aoa > alphaL
    m = (CL_H-CL_L)/(alphaH - alphaL);
    b = CL_L - m*alphaL;
    coeffLift = m*aoa + b;
    coeffDrag = profileDrag + coeffLift^2/(pi*oswaldE*AR);
    lastLinearCD = coeffDrag;
elseif aoa > alphaH && aoa <= (alphaH+8)
    m = (CL_H- coeffLDEnd)/(alphaH - (8+alphaH));
    coeffLift = m*aoa + CL_H - m*alphaH;
    coeffDrag = lastLinearCD + (lastLinearCD-coeffDHEnd)/(alphaH - (alphaH+8))*(aoa-10);
else
    coeffLift = cos(2.0*((aoa/180*pi)-45/180*pi));
    coeffDrag = sin(2*((aoa-45)/180*pi))+1.0;
end
LiftDrag(1) = coeffLift;
LiftDrag(2) = coeffDrag;
end

function launchForce = getLaunchForce(t)
g = 9.98;
tLimit = 1;
if t < tLimit
    launchForce = 20*g - 20*g*t;
else
    launchForce = 0;
end
end

I have no idea what the expected results are, but here is what this code produced for me.

[Megan Frisbey]

Thank you.

"return a vector of the 12 values at that time step"

Okay yes, that is what I was originally trying to do (sort of) by returning a vector of 12 symbolic ODE's unique to that time step (really its only verying a constant for a force term at that time step)

function dydt = osc_1BodyKane_2Verification011823(time,Y)

(same stuff here)

%Y = [q4, Dq4, q1, Dq1, q6, Dq6, q3, Dq3, q5, Dq5, q2, Dq2]
[dydt,Y] = odeToVectorField(Kanes1 == 0,Kanes2 == 0, Kanes3 == 0, Kanes4 == 0, Kanes5 == 0, Kanes6 == 0);
end

Which would complain about

Error using superiorfloat

Inputs must be floats, namely single or double

Coincidentally I get the same error when I try to implement your suggested changes above :( I don't have all my functions in the same file however. Could this be a scope issue? I'm not familiar enough with how matlab does scoping and global variables. I prefer not to have all the code in the same file. Also when I hand code the output from odeToVectorField from say the first time step as below it will work okay sending the dydt vector back (see below, but it's ugly to look at and only a sanity check)

% dydt(1) = y(2);
% dydt(2) = ((cos(y(9)).*sin(y(9)).*y(2).*y(10).*2.3005e+4+cos(y(1)).*cos(y(9)).^2.*sin(y(1)).*y(6).^2.*1.1124e+4-cos(y(1)).*sin(y(1)).*sin(y(9)).^2.*y(6).^2.*1.1881e+4+cos(y(1)).*cos(y(9)).^2.*y(6).*y(10).*2.4192e+4+cos(y(1)).*sin(y(9)).^2.*y(6).*y(10).*4.7197e+4-cos(y(9)).*sin(y(1)).*sin(y(9)).*y(2).*y(6).*2.3005e+4).*(-2.831577755125156e-5))./(cos(y(9)).^2+sin(y(9)).^2);
% dydt(3) = y(4);
% dydt(4) = cos(y(9)).^2.*y(2).^2.*y(3)+sin(y(9)).^2.*y(2).^2.*y(3)-cos(y(9)).*y(2).*y(8).*2.0-sin(y(9)).*y(2).*y(12).*2.0+cos(y(1)).*cos(y(9)).*y(6).*y(12).*2.0-cos(y(1)).*sin(y(9)).*y(6).*y(8).*2.0+cos(y(1)).^2.*cos(y(9)).^2.*y(3).*y(6).^2+cos(y(1)).^2.*sin(y(9)).^2.*y(3).*y(6).^2-cos(y(9)).*y(2).*y(10).*y(11).*(2.15e+2./3.24e+2)+sin(y(9)).*y(2).*y(7).*y(10).*(4.3e+2./3.27e+2)+cos(y(1)).*cos(y(9)).*sin(y(1)).*y(6).^2.*y(7).*(4.3e+2./3.27e+2)+cos(y(1)).*sin(y(1)).*sin(y(9)).*y(6).^2.*y(11).*(2.15e+2./3.24e+2)-cos(y(1)).*cos(y(9)).*y(6).*y(7).*y(10).*(4.3e+2./3.27e+2)+cos(y(9)).*sin(y(1)).*y(2).*y(6).*y(11).*(2.15e+2./3.24e+2)-cos(y(1)).*sin(y(9)).*y(6).*y(10).*y(11).*(2.15e+2./3.24e+2)-sin(y(1)).*sin(y(9)).*y(2).*y(6).*y(7).*(4.3e+2./3.27e+2);
% dydt(5) = y(6);
% dydt(6) = (cos(y(9)).^2.*y(2).*y(10).*4.7197e+4+sin(y(9)).^2.*y(2).*y(10).*2.4192e+4+cos(y(9)).^2.*sin(y(1)).*y(2).*y(6).*2.3435e+4+sin(y(1)).*sin(y(9)).^2.*y(2).*y(6).*4.644e+4-cos(y(1)).*cos(y(9)).*sin(y(1)).*sin(y(9)).*y(6).^2.*2.3005e+4+cos(y(1)).*cos(y(9)).*sin(y(9)).*y(6).*y(10).*2.3005e+4)./(cos(y(1)).*sin(y(9)).^2.*3.5316e+4+cos(y(1)).*cos(y(9)).^2.*3.5316e+4);
% dydt(7) = y(8);
% dydt(8) = y(7).*y(10).^2-y(10).*y(12).*2.0+cos(y(9)).^2.*y(2).^2.*y(7)+sin(y(1)).^2.*y(6).^2.*y(7)+cos(y(9)).*y(2).*y(4).*2.0+sin(y(1)).*y(6).*y(12).*2.0+cos(y(1)).*sin(y(9)).*y(4).*y(6).*2.0+cos(y(1)).^2.*sin(y(9)).^2.*y(6).^2.*y(7)+cos(y(9)).*sin(y(9)).*y(2).^2.*y(11).*(2.24e+2./5.45e+2)-sin(y(1)).*y(6).*y(7).*y(10).*2.0+sin(y(9)).*y(2).*y(3).*y(10).*(2.24e+2./3.27e+2)-cos(y(1)).*cos(y(9)).^2.*y(2).*y(6).*y(11).*(2.24e+2./5.45e+2)+cos(y(1)).*sin(y(9)).^2.*y(2).*y(6).*y(11).*(2.24e+2./5.45e+2)+cos(y(1)).*cos(y(9)).*sin(y(1)).*y(3).*y(6).^2.*(2.24e+2./3.27e+2)-cos(y(1)).*cos(y(9)).*y(3).*y(6).*y(10).*(2.24e+2./3.27e+2)-sin(y(1)).*sin(y(9)).*y(2).*y(3).*y(6).*(2.24e+2./3.27e+2)-cos(y(1)).^2.*cos(y(9)).*sin(y(9)).*y(6).^2.*y(11).*(2.24e+2./5.45e+2)+cos(y(1)).*cos(y(9)).*sin(y(9)).*y(2).*y(6).*y(7).*2.0;
% dydt(9) = y(10);
% dydt(10) = (cos(y(1)).*cos(y(9)).*sin(y(9)).^3.*y(2).^2.*1.04004e+5+cos(y(1)).*cos(y(9)).^3.*sin(y(9)).*y(2).^2.*1.04004e+5+cos(y(9)).^2.*sin(y(1)).*y(2).*y(10).*2.35985e+5+sin(y(1)).*sin(y(9)).^2.*y(2).*y(10).*1.2096e+5-cos(y(1)).^3.*cos(y(9)).*sin(y(9)).^3.*y(6).^2.*1.04004e+5-cos(y(1)).^3.*cos(y(9)).^3.*sin(y(9)).*y(6).^2.*1.04004e+5+cos(y(1)).^2.*cos(y(9)).^2.*y(2).*y(6).*1.7658e+5-cos(y(1)).^2.*cos(y(9)).^4.*y(2).*y(6).*1.04004e+5+cos(y(1)).^2.*sin(y(9)).^2.*y(2).*y(6).*1.7658e+5+cos(y(9)).^2.*sin(y(1)).^2.*y(2).*y(6).*1.17175e+5+cos(y(1)).^2.*sin(y(9)).^4.*y(2).*y(6).*1.04004e+5+sin(y(1)).^2.*sin(y(9)).^2.*y(2).*y(6).*2.322e+5-cos(y(1)).*cos(y(9)).*sin(y(1)).^2.*sin(y(9)).*y(6).^2.*1.15025e+5+cos(y(1)).*cos(y(9)).*sin(y(1)).*sin(y(9)).*y(6).*y(10).*1.15025e+5)./(cos(y(1)).*sin(y(9)).^2.*1.7658e+5+cos(y(1)).*cos(y(9)).^2.*1.7658e+5);
% dydt(11) = y(12);
% dydt(12) = y(10).^2.*y(11)+y(8).*y(10).*2.0+sin(y(1)).^2.*y(6).^2.*y(11)+sin(y(9)).^2.*y(2).^2.*y(11)-sin(y(1)).*y(6).*y(8).*2.0+sin(y(9)).*y(2).*y(4).*2.0-cos(y(1)).*cos(y(9)).*y(4).*y(6).*2.0+cos(y(1)).^2.*cos(y(9)).^2.*y(6).^2.*y(11)+cos(y(9)).*sin(y(9)).*y(2).^2.*y(7).*(8.66e+2./5.45e+2)-cos(y(9)).*y(2).*y(3).*y(10).*(4.33e+2./3.24e+2)-sin(y(1)).*y(6).*y(10).*y(11).*2.0-cos(y(1)).*cos(y(9)).^2.*y(2).*y(6).*y(7).*(8.66e+2./5.45e+2)+cos(y(1)).*sin(y(9)).^2.*y(2).*y(6).*y(7).*(8.66e+2./5.45e+2)+cos(y(1)).*sin(y(1)).*sin(y(9)).*y(3).*y(6).^2.*(4.33e+2./3.24e+2)+cos(y(9)).*sin(y(1)).*y(2).*y(3).*y(6).*(4.33e+2./3.24e+2)-cos(y(1)).*sin(y(9)).*y(3).*y(6).*y(10).*(4.33e+2./3.24e+2)-cos(y(1)).^2.*cos(y(9)).*sin(y(9)).*y(6).^2.*y(7).*(8.66e+2./5.45e+2)-cos(y(1)).*cos(y(9)).*sin(y(9)).*y(2).*y(6).*y(11).*2.0;

and I'm at a loss to determine the difference between what I hand code in and what odeToVectorField is returning in regards to data type and why one would be okay and the other wouldn't? Is it because I'm using syms q1(t)... etc etc

Also just to note I only uncommented some of the output for debugging purposes! That would defintely help with speed up but an hour run time is pretty terrible compared to what I'm seeing with other methods.

Once again thank you for a thoughtful reply!

[Cris LaPierre]

It's not a variable scope related issue. The way I have written them, the behavior should be the same as if each function were a separate file.

BTW, tt took 25 minutes to run with all outputs suppressed.

[Megan Frisbey]

Okay I found the difference between when I tried to implement your change and then just copy pasting your corrected code. Can you help me understand the difference between these?

dydt = odeToVectorField(Kanes1 == 0,Kanes2 == 0, Kanes3 == 0, Kanes4 == 0, Kanes5 == 0, Kanes6 == 0);
odefun = matlabFunction(dydt,'vars', {t,'Y'});
Ynew = odefun(time,Y); %################  Add this line to evaluate odefun at current time step

I tried this

[dydt,Y] = odeToVectorField(Kanes1 == 0,Kanes2 == 0, Kanes3 == 0, Kanes4 == 0, Kanes5 == 0, Kanes6 == 0);
odefun = matlabFunction(dydt,'vars', {t,'Y'});
Ynew = odefun(time,Y); %################  Add this line to evaluate odefun at current time step

why does the latter give the error

Error using superiorfloat

Inputs must be floats, namely single or double

and the one above "works" Y is a vector of the subs?

25 minutes is improved over an hour but I have other methods that run in about 10 seconds. Is it the symbolic manipulation? the odeToVectorField seems to eat up quite a bit of time.

Once again thank you for your very helpful responses

[Cris LaPierre]

In the second code snippet, odeToVectorField is returning a second output, Y. This is the same variable name as your 2nd input to the osc_1BodyKane_2Verification011823 function, and therefore overwrites it. When you go to the next line of code, Y is no longer a vector of numbers, but is instead now a vector of symbolic substitutions. This is producing the error.

So there is nothing wrong with the syntax of getting a second output from odeToVectorField. You just need to use a different variable name for it.

[dydt,Ysubs] = odeToVectorField(Kanes1 == 0,Kanes2 == 0, Kanes3 == 0, Kanes4 == 0, Kanes5 == 0, Kanes6 == 0);

But since you don't need or use that output, you can also drop it entirely, as I did.