# 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.

5 views (last 30 days)
Megan Frisbey on 24 Jan 2023
Edited: Cris LaPierre on 25 Jan 2023
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!
Megan Frisbey on 25 Jan 2023
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.

Cris LaPierre on 24 Jan 2023
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
##### 2 CommentsShowHide 1 older comment
Megan Frisbey on 25 Jan 2023
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

Cris LaPierre on 25 Jan 2023
Edited: Cris LaPierre on 25 Jan 2023
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.
Cris LaPierre on 25 Jan 2023
Edited: Cris LaPierre on 25 Jan 2023
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.

### Categories

Find more on LMI Solvers in Help Center and File Exchange

R2022b

### Community Treasure Hunt

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

Start Hunting!