Problems Implementing ddesd with a system of equations.
5 次查看(过去 30 天)
显示 更早的评论
I am attempting to model a system of equations to capture a single rigid body dynamics problem. I want to use functions for the forces and moments acting on the body thus creating a state dependent function. First question - this does make it a DDE problem correct? Next I have implemented a simple rudemetary model using the examples given online but I am missing/not understanding something when it comes to the history function. I took this as a way to define the initial conditions, but after seeing the results I believe this is not correct or I have not defined it properly? To back out what was going on I dusted off my ODE solver using ode45 with the same equations and such and determined that the ddesd solver is using the history function and forcing those values for my delay values -- Why?!? It doesn't matter how I change my delay functions it has no effect on the solution and appears to be only using the value for the first position of the vector I define in the histroy function. Obvously I'm missing something, but what?
I appreciate the help!
below is the dde system I've tried to implement.
t0 = 0
tf = 10
t = linspace(t0,tf,200)
sol = ddesd(@ddefun,@delay,@history,t)
t = sol.x
y = sol.y
%I've ommited all the plotting but can add this back in if needed
function d = delay(t,y)
d = [delayForce(t,y) delayMoment(t,y)]
end
function dF = delayForce(t,y)
g = 32.174 %ft/s^2
mB = 50; %lbs
rho = 0.002377 %SL english
% FLaunch = getLaunchForce(t);
% FAero = getAeroForce(y,rho)
% FWeight = mB %mB*g
% FThrust = -FAero(1)*20
% Fbx = FLaunch + FAero(1) + FThrust;
% Fby = FAero(2);
% Fbz = FAero(3) + -FWeight;
Fbx = 0; %getLaunchForce(t)
Fby = 0;
Fbz = 0;
dF = [Fbx Fby Fbz];
end
function dM = delayMoment(t,y)
Lbx = 0
Mbx = 0
Nbx = 0
dM = [Lbx Mbx Nbx]
end
function v = history(t)
% U = y(1); V = y(2); W = y(3); x = y(7); y = y(8); z = y(9);
% P = y(4); Q = y(5); R = y(6); phi = y(10); theta = y(11); psi = y(12);
g = 32.174 %ft/s^2
v = transpose([10 0 0 0 0 0 0 0 0 0 0 0]);
end
function dydt = ddefun(t,y,Z)
dydt = zeros(12,1);
mB = 50/32.174; %lbs
mPlate = 6;
rPlate = 6;
hPlate = 1;
IBxx = 1/12*mB*(3*rPlate*rPlate + hPlate*hPlate);
IByy = 1/12*mPlate*(3*rPlate*rPlate + hPlate*hPlate);
IBzz = 1/2*mPlate*rPlate*rPlate;
% u1 = y(1); u2 = y(2); u3 = y(3); q1 = y(7); q2 = y(8); q3 = y(9);
% u4 = y(4); u5 = y(5); u6 = y(6); q4 = y(10); q5 = y(11); q6 = y(12);
dydt(1) = Z(1,1)*(1/mB) + y(2)*y(6) - y(3)*y(5);
dydt(2) = Z(1,2)*(1/mB) + y(3)*y(4) - y(1)*y(6);
dydt(3) = Z(1,3)*(1/mB) + y(1)*y(5) - y(2)*y(4);
dydt(4) = (Z(1,4) + (IByy - IBzz)*y(5)*y(6))*(1/IBxx);
dydt(5) = (Z(1,5) + (IBzz - IBxx)*y(4)*y(6))*(1/IByy);
dydt(6) = (Z(1,6) + (IBxx - IByy)*y(4)*y(5))*(1/IBzz);
dydt(7) = y(1) + y(6)*y(8) - y(5)*y(9);
dydt(8) = y(2) + y(4)*y(9) - y(6)*y(7);
dydt(9) = y(3) + y(5)*y(7) - y(4)*y(8);
dydt(10) = y(4) + y(5)*sin(y(10))*tan(y(11)) +y(6)*cos(y(10))*tan(y(11));
dydt(11) = y(5)*cos(y(10)) -y(6)*sin(y(10));
dydt(12) = y(6)*(cos(y(10))/cos(y(11)))+y(5)*(sin(y(10))/cos(y(11)));
end
Below here is the code where I use ode45 and having tuned the force functions so that it matches the incorrect output from ddesd.
function [T,Y] = call_osc_1BodyNE()
tspan = [0 10];
initVec = [10 0 0 0 0 0 0 0 0 0 0 0]
y1_0 = 2;
y2_0 = 0;
[T,Y] = ode45(@osc_1BodyNE_alt,tspan,initVec);
%I've ommited all the plotting but can add this back in if needed
function dydt = osc_1BodyNE_alt(t,y)
dydt = zeros(12,1);
g = 32.174;
weight = 50;
rho = 0.002377;
massB = weight/g;
mB = massB %Generic
mPlate = 6;
rPlate = 6;
hPlate = 1;
%Calculate Forces
% FLaunch = getLaunchForce(t);
% FAero = getAeroForce(y,rho)
% FWeight = massB*g
% FThrust = FAero(1)
%
% F_Applied(1) = FLaunch + FAero(1)+FThrust;
% F_Applied(2) = FAero(2);
% F_Applied(3) = FAero(3) + -FWeight;
% F_Applied(1) = getLaunchForce(t);
F_Applied(1) = 10;
F_Applied(2) = 10;
F_Applied(3) = 10;
% + FAero(3)*1.2; %
% Lift = FAero(3)
%Mass Moments
Ib = zeros(3,3);
Ib(1,1) = 1/12*massB*(3*rPlate*rPlate + hPlate*hPlate);
Ib(2,2) = 1/12*mPlate*(3*rPlate*rPlate + hPlate*hPlate);
Ib(3,3) = 1/2*mPlate*rPlate*rPlate;
IBxx = Ib(1,1);
IByy = Ib(2,2);
IBzz = Ib(3,3);
%Calculate Moments
M_Applied(1) = 10;
M_Applied(2) = 10;
M_Applied(3) = 10;
%body
% U = y(1); V = y(2); W = y(3); x = y(7); y = y(8); z = y(9);
% P = y(4); Q = y(5); R = y(6); phi = y(10); theta = y(11); psi = y(12);
dydt(1) = F_Applied(1)*(1/mB) + y(2)*y(6) - y(3)*y(5);
dydt(2) = F_Applied(2)*(1/mB) + y(3)*y(4) - y(1)*y(6);
dydt(3) = F_Applied(3)*(1/mB) + y(1)*y(5) - y(2)*y(4);
dydt(4) = (M_Applied(1) + (IByy - IBzz)*y(5)*y(6))*(1/IBxx);
dydt(5) = (M_Applied(2) + (IBzz - IBxx)*y(4)*y(6))*(1/IByy);
dydt(6) = (M_Applied(3) + (IBxx - IByy)*y(4)*y(5))*(1/IBzz);
dydt(7) = y(1) + y(6)*y(8) - y(5)*y(9);
dydt(8) = y(2) + y(4)*y(9) - y(6)*y(7);
dydt(9) = y(3) + y(5)*y(7) - y(4)*y(8);
dydt(10) = y(4) + y(5)*sin(y(10))*tan(y(11)) +y(6)*cos(y(10))*tan(y(11));
dydt(11) = y(5)*cos(y(10)) -y(6)*sin(y(10));
dydt(12) = y(6)*(cos(y(10))/cos(y(11)))+y(5)*(sin(y(10))/cos(y(11)));
end
回答(1 个)
Sulaymon Eshkabilov
2023-6-28
Why you are not using the time span with a predefined time step, in call_osc_1BodyNE(), that would give you consistent time steps in your found solutions from ode45 e.g.:
function [T,Y] = call_osc_1BodyNE()
t0 = 0;
tf = 10;
tspan = linspace(t0,tf,200);
initVec = [10 0 0 0 0 0 0 0 0 0 0 0];
y1_0 = 2;
y2_0 = 0;
[T,Y] = ode45(@osc_1BodyNE_alt,tspan,initVec);
%I've ommited all the plotting but can add this back in if needed
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!