Clear Filters
Clear Filters

Problems Implementing ddesd with a system of equations.

10 views (last 30 days)
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 Comments
Megan Frisbey
Megan Frisbey on 29 Jun 2023
In my question I specifically referred to using ode45 which is only dependent on the previous time step
"ode45 is based on an explicit Runge-Kutta (4,5) formula, the Dormand-Prince pair. It is a single-step solver – in computing y(tn), it needs only the solution at the immediately preceding time point, y(tn-1) [1], [2]."
Torsten
Torsten on 29 Jun 2023
Edited: Torsten on 29 Jun 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:

Sign in to comment.

Answers (1)

Sulaymon Eshkabilov
Sulaymon Eshkabilov on 28 Jun 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 Comments
Megan Frisbey
Megan Frisbey on 29 Jun 2023
No good reason I suppose other than it isn't really necessary? Is there a compelling reason to do this using ode45? The solution doesn't change and ode45 seems to be doing fine managing the time step...
Megan Frisbey
Megan Frisbey on 29 Jun 2023
This is not an answer. Did you "answer" this question just so it would stay around in matlab answers archieves?

Sign in to comment.

Tags

Products


Release

R2022b

Community Treasure Hunt

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

Start Hunting!