Problems Implementing ddesd with a system of equations.
5 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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
27 comentarios
Torsten
el 29 de Jun. de 2023
Editada: Torsten
el 29 de Jun. de 2023
That's true. But it will call your ode function not only for t(n-1) and y(t(n-1)) to compute y(t(n)), but also for some "intermediate" points. You can see this here for the classical Runge-Kutta method - it's similarily done in ode45:
Respuestas (1)
Sulaymon Eshkabilov
el 28 de Jun. de 2023
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
2 comentarios
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!