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.

3 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!
5 CommentsShow 3 older commentsHide 3 older comments
Megan Frisbey on 25 Jan 2023
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 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.

Sign in to comment.

Answers (2)

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 CommentsShow NoneHide None
Cris LaPierre on 24 Jan 2023
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 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

Sign in to comment.

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.
4 CommentsShow 2 older commentsHide 2 older comments
Megan Frisbey on 25 Jan 2023
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 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.

Sign in to comment.

Categories

Find more on Ordinary Differential Equations 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!