Problems Implementing ddesd with a system of equations.

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

[Torsten]

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?

"State-dependent" does not mean "work with a time delay". So why do you think you need to incorporate a temporal delay in your equations ?

[Megan Frisbey]

The forces and moments are a function of the orientation and velocity (part of the y solution) from the previous step. So I assumed this made it a DDE not simply an ODE?

[Torsten]

And how long is the previous step ? 1/10 s, 1/100 s, 1/1000 s ? In solving ODEs, something like "the previous step" does not exist.

[Megan Frisbey]

Previous time step? ode45 uses a runge-kutta scheme which has different steps or half steps or some deriviative of a step depending on the order, so that response seems misleading. Regardless my question is I have a system of equations one of the inputs for several of those equations is not a constant or linear function, but rather a function of the previous state. I asked previously about ode45 and solving this set of what I thought at the time were ode's and from this forum was told that it was probably a dde problem that I should consider different options? (https://www.mathworks.com/matlabcentral/answers/1900090-implementation-issues-in-function-with-symbolic-ode-using-ode45-name-returns-a-vector-of-length-1-bu#comment_2581110)

Having never worked with dde's before I was hoping to have confirmation on this. And then of course some insite into why the history function is behaving the way it is. Thanks!

[Torsten]

Z(1,1) means y(1), evaluated at the first delay.

Z(1,2) means y(1), evaluated at the second delay.

...

Z(1,6) means y(1), evaluated at the sixth delay.

Your delays are all defined as 0.

Thus Z(1,1) = Z(1,2) = ... = Z(1,6) = y(1), evaluated at t = 0 for all times t.

Does that help ?

[Megan Frisbey]

No, - I defined my delays as zero for debugging purposes. This allowed me to see that the value of 10 in the history vector was being input for all of the delay values as well. If I set the delay terms back to their functions (commented out above), the history vector acts as a mulitplier. I don't understand this behavior. From the documentation the history vector is supposed to be the state at time less than initial time. How/why is this being factored in? Should I be setting my Z terms in the equations to Z(t,1) Z(t,2) etc?

Also to confirm this is or is not a dde system?

[Torsten]

The delays are set as

t - delaylength

not as

delaylength

as you do. I think you didn't consider this, did you ?

In the example provided under

https://uk.mathworks.com/help/matlab/ref/ddesd.html

e.g.

function d = ddex1delays(t,y)
%DDEX1DELAYS  Delays for using with DDEX1DE.
d = [ t - 1
    t - 0.2];
end

[Megan Frisbey]

Edit to add, yes that maybe be helpful some. But I'm also a little confused. Isn't the first term in Z the timestep? y is a vector here I assumed the whole vector was passed as I need mulitple components of y to evaluate the delay functions for each equation.

[Torsten]

Example:

If you specify

d(1) = t-pi

, then

Z(1,1) = y1(t-pi)
Z(2,1) = y2(t-pi)
...
Z(n,1) = yn(t-pi)

where y1(t-pi),...,yn(t-pi) are your solution variables, all evaluated at t-pi (and t is the time instant passed to "ddefun").

[Megan Frisbey]

Yes I've seen that example. I didn't find it very useful.

I am missing (clearly lol) a major concept here. Does the delay function not simply return a value for the term Z(t,y) set in the equations for the ddefun function? In that example given it appear that instead it is returning the y values at some other time. I have used this example to try and follow instead since I am state dependent not time dependent.

https://www.mathworks.com/help/matlab/math/state-dependent-delay-problem.html

function d = dely(t,y)
d = exp(1 - y(2));
end

where the function dely returns the value d, not y evaluated at a different time.

[Megan Frisbey]

Also just to note there is a delay in my responses and your responses. The last comment I made was in reference to this comment. https://www.mathworks.com/matlabcentral/answers/1988933-problems-implementing-ddesd-with-a-system-of-equations#comment_2798398

I think some of the confusion coming in is around the type of delay. I have state dependent delays, where my delay functions depend on my solution from the previous time solution. So once again the first question - Is this truely a DDE system? Secondly how is the history term supposed to behave. (I have seen the documention here https://www.mathworks.com/help/matlab/math/state-dependent-delay-problem.html but this is not what I'm observing in my solution)

[Torsten]

Does the delay function not simply return a value for the term Z(t,y) set in the equations for the ddefun function?

No. The delay function returns the time argument at which the solution variables y1,...,yn are to be evaluated in "ddefun".

In that example given it appear that instead it is returning the y values at some other time.

Yes, this is exactly the principle of delay differential equations.

Did you look at the example here:

https://uk.mathworks.com/help/matlab/math/state-dependent-delay-problem.html

?

This is a state-dependent delay problem. The function y2 in "ddefun" is to be evaluated at time exp(1-y2), not at time t. That's why d(1) is set to exp(1-y(2)).

[Megan Frisbey]

That is the exact example I just linked and referenced? So I completely misunderstood what the delay function was doing and really the over all problem (facepalm). This leads me back around to either I need to reformulate my equations to incorperate a constant time delay (to get the solution from the previous step) OR this really isn't a DDE problem...any suggestions? Am I tracking now? Thanks!

[Torsten]

I don't know what your equations are in a mathematical notation and what you try to model. Thus I'm not sure what you really mean by "solution from the previous step" (because "the previous step" is not fixed in size and thus arbitrary) and whether this is the correct way to handle a certain term appearing in your system of equations.

In the equations you gave the link to

https://www.mathworks.com/matlabcentral/answers/1900090-implementation-issues-in-function-with-symbolic-ode-using-ode45-name-returns-a-vector-of-length-1-bu#comment_2581110

I cannot see a term that uses the solution of a previous time step.

[Megan Frisbey]

Maybe I need to reframe the problem for you to better communicate what I mean. This is an IVP the previous step I'm refering to is the y solution of dydt that was generate for (time-deltat) (or in the case of the first solution at t0 the initial conditions) The size of the time step for this is arbitrary (assuming stability) in a time marching scheme each slice of time is what I'm referring to as "steps" perhaps iterations, counts, solution at time t? If I'm solving at time t I need the solution at time (t-dt) to serve as my initial values and my force and moments are dependent on this also.

From the link above you will find the equations that use the previous time step solution in this funciton, y being the solution found at the previous iteration.

function AeroForce = getAeroForce(y,rho)

[Torsten]

Although I doubt that your mathematical ansatz of using "the solution at the previous time step" is correct, you have two options: either use dde23 with the (fixed) delay defined as the deltat that you think is adequate or use a fixed time step solver where you have access to the solution at previous times.

[Megan Frisbey]

Okay can you explain what you think is wrong about this approach? I've derived the basic newton-euler 6dof equations of motion for a single rigid body. This body is airborne and the forces and moments are derived from the aerodynamic effects that are dependent on the bodies velocity and orientation (part of the solution of the set of equations)

[Torsten]

I don't know about the physics. But I've never seen a mathematical model where the equations refer to the values of the solution variables of the preceding time step, especially if this time step is not quantified. You still didn't explain why you can't work with the instantaneous values of the solution variables.

[Megan Frisbey]

Obtaining accurate solutions for fluid-structure interaction problems is very complex. The equations that rule CFD is the Navier Stokes equations - PDE's, and also for structures more PDE's. I am attempting a much more simplified less refined method of doing this by precomputing the values of coefficient of drag and lift on a body at different flow field states, and by only considering the rigid body dyanimcs (no bending etc). For a very crude model I can reduce this down to some peice wise function. For a better model it will be table look up based off of conditions. I could attempt to map the solutions surface and formulate into some form of an equation but this I believe is more trouble than it is worth. If the time step is small enough for this problem using the previous state values to determine the current force and moment values would be "good enough"

As to the time step not being quantified I'm not sure what you mean. Plenty of numerical schemes don't quanitfy a time step outside of stability bounds and accuracy of solution...What are you meaning by this?

[Torsten]

If the time step is small enough for this problem using the previous state values to determine the current force and moment values would be "good enough"

And why can't you use the actual state values to determine the current force and moment values ? Are they computed in different programs ?

[Megan Frisbey]

The actual state values are what I am solving for???

[Megan Frisbey]

And the forumulation of those variables in the force and moment calculations is too complicated to just wrap up in my set of equations.

[Torsten]

The actual state values are what I am solving for???

Yes, and they make up the y-vector that is handed to you in your ODE function by every MATLAB integrator so that you can use them to calculate all derived quantities (such as force and moment).

[Megan Frisbey]

Okay so let me translate what you are saying -

1) you are handed the y-vector at (t-dt) previous time step solution so why not also use this?

Oh my goodness - this is so circular. It's the very first question I asked.....is this a DDE problem? 20 something comments in -what you are suggesting is that it is not because the functions don't have a dependency on a solution other than the previous time step.

[Torsten]

It's not necessarily true that you are handed the solution at time t-dt. For the Explicit Euler method with fixed time stepping, this will be the case.

What you get for t and y in the function depends on the discretization scheme for the solution method used to solve your differential equations (Euler explicit, Euler implicit, Runge-Kutta,...). But t and y are the values with which you have to evaluate the function f of your system dy/dt = f(t,y) and from which you have to deduce derived quantities.

[Megan Frisbey]

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]

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:

https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

Answer 1

[Sulaymon Eshkabilov]

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

[Megan Frisbey]

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]

This is not an answer. Did you "answer" this question just so it would stay around in matlab answers archieves?