Update of Variables during the use of ODE113

Good afertnoon,

I am trying to implement a piece of code to solve Multibody Dynamics Simualtions. I pretend to use the the solver ode113, but I am having some trouble on how to manage the variables.

My question is the following, at each time step/integration of the solver, I have to update my values of position and velocity, this update should be done inside the odefun after I give yd or at the start of the function? Is this update of the variables going to affect the accuracy of the odesolver? This is, will it update with the predictor steps or only after the tolerances are met?

This is my current code for the odefun (which is not 100% finished due to my doubts):

%% Set Up of the variables needed to construct the Matrix - Dynamic Modified Jacobian
% Pre Allocation of the q0 vector, for the calc of AGL
    q0 = CreateAuxiliaryBodyStructure(NBodies,Bodies);    
% Calculation of the Matrixes A, G and L
    Bodies = DynCalcAGL(q0,NBodies,Bodies);
% Calculus of the dofs
    coord = 7;
    [debugdata] = SystemDofCalc(NBodies,Joints,debugdata,[],coord);
% Calculus of the G and L matrices derivatives
    [Bodies] = DynIAGdLdcalc(NBodies,Bodies);
% Change w field to allow Kinematic Functions to Run on the RHS Accel (Can
% be Optimized)

    for i = 1:NBodies
       Bodies(i).wl = Bodies(i).w; 
    end
    
%% Assembly of the Invariable Matrices in the Augmented Lagrangian Formula.
    % The process to obtain Position Matrix will be separated from the
    % Jacobian to avoid funcount conflitcs
    
    %Position Matrix
    Flags.Position = 1;
    Flags.Jacobian = 0;
    Flags.Velocity = 0;
    Flags.Acceleration = 0;
    Flags.Dynamic = 0;
    Flags.AccelDyn = 0;
    
    [fun,~,~,~] = PJmatrixfunct(Flags,Bodies,NBodies,Joints,debugdata,driverfunctions,coord);
    
    %Jacobian Matrix
    Flags.Position = 0;
    Flags.Jacobian = 1;
    Flags.Velocity = 0;
    Flags.Acceleration = 0;
    Flags.Dynamic = 0;
    Flags.AccelDyn = 0;
    
    [~,Jacobian,~,~] = PJmatrixfunct(Flags,Bodies,NBodies,Joints,debugdata,driverfunctions,coord);
    
%% Mass Matrix Assembly for Euler Parameters
    massmatrix = zeros(7*NBodies,7*NBodies); %Pre-Allocation
    for i = 1:NBodies
        Mass = Bodies(i).Mass;
        B = Bodies(i).L;
        Inertia = Bodies(i).Inertia;
        Irat = 4*B'*diag(Inertia)*B; %Nikra Article on the use of EP for 3D Dynamics
        i1 = 7*(i-1)+1;
        massmatrix(i1:i1+2,i1:i1+2) = Mass * eye(3);
        massmatrix(i1+3:i1+6,i1+3:i1+6) = Irat; 
    end
    
%% Definition of the Velocity Vector from 6 coordinates to 7 coordinates
    for i = 1:NBodies
        i1 = 7*(i-1)+1;
        qd(i1:i1+2,1) = Bodies(i).rd;
        pd = 0.5*Bodies(i).L'*Bodies(i).w;
        pd = Impose_Column(pd);
        qd(i1+3:i1+6,1) = pd;
    end
    
%% Function Responsible for the Force Vectors
    [vetorg] = ForceCalculus(Forces,NBodies,Bodies,dynfunc,Grav,UnitsSystem,time,ForceFunction);

%% Solving First Iteration to start the Augmented Process
    [dim1,dim2] = size(massmatrix);
    wogmass = massmatrix(8:dim1,8:dim2);
    [dim3,~] = size(vetorg);
    wogvetor = vetorg(8:dim3,1);
    qddwog = pinv(wogmass)*wogvetor;
    [dim4,~] = size(qddwog);
    qdd0(8:dim4+7,1) = qddwog; 

%% Solving the Augmented Lagrangian Formula Iterative Formula
    alf = 'alfon';
    beta = 1;
    qddi = qdd0;
    lagit = 1; % Augmented Lagrangian Formula Iteration Numver
    % Flags to retrieve gamma
    Flags.Position = 0;
    Flags.Jacobian = 0;
    Flags.Velocity = 1;
    Flags.Acceleration = 1;
    Flags.Dynamic = 0;
    Flags.AccelDyn = 0;
    while beta > 1e-6 % Second Condition to avoid infinite loops
        if lagit > 1
           qddi = qddi1; 
        end
        [~,~,niu,gamma] = PJmatrixfunct(Flags,Bodies,NBodies,Joints,debugdata,driverfunctions,coord);
        lhslag = massmatrix + Jacobian'*alpha*Jacobian;
        rhslag = massmatrix*qddi + Jacobian'*alpha*(gamma - 2*omega*mu*(Jacobian*qd - niu) - omega^2*fun);
        qddi1 = pinv(lhslag)*rhslag;
        beta = alpha*(Jacobian*qddi1 - gamma + 2*omega*mu*(Jacobian*qd - niu) + omega^2*fun);
        [Bodies] = UpdateAccelerations(qddi1,NBodies,Bodies,SimType,alf);
        lagit = lagit + 1; %Iteration Counter
    end

    yd = [qd;qddi1];
    
    

[Torsten]

I don't understand exactly what you are asking for.

The functions in which you provide the derivatives and/or the mass matrix get t and y(j) as inputs.

You have to calculate the derivatives dy(j) and/or the mass matrix elements A(i,j) that follow from these input data.

[Tiago Carvalho]

Sorry if I wasn't explicit enough. I know that I have to input the derivatives (velocity and acceleration) to allow the solver to integrate and give me the (position and velocity). However, at each time step, my data structures which contain the information needed to calculate those derivatives need to be updated with the previous iteration output (position and velocity) and my difficulty is knowing when to update these variables on my odefun.

[Torsten]

Usually, you can only update your data structures with the values of the current iteration because you don't have access to the previous iterations's values. It's not true for MATLAB's ODE solvers that call i of your function is for t(i) and call (i+1) is for time t(i+1) where t is the array of output times you defined in tspan in the call to ODE113. MATLAB's time integration is adaptive - so there might be many intermediate calls to your function.

The only thing you can do is to program a function OutputFcn that is called after each successful time step, update your data and work with these data within the time period to the next successful time step. But this can cause inexact result and even convergence issues, of course.

[Tiago Carvalho]

But for example, it is not possible to call the results y from the previous integration step to the odefun?

[Tiago Carvalho]

I only have photos of a code that tries to accomplish the same:

Script:

odefun:

[Tiago Carvalho]

From what I could gather from this code, the inputs of the odefun are calling t,y that are also the outputs of the integration in order to update, or not?

[Torsten]

No, "odefun" only works with "trial" values for y for the actual time t and does not provide values after the last successful time step.

As I said, you can update your structure data in "OutputFcn" with the results of the actual time t and use these data in the time interval from t to the next successful time step t+dt.

Look at "OutputFcn" under

https://de.mathworks.com/help/matlab/ref/odeset.html

for more details.

[Tiago Carvalho]

Ok ok, I think I got it thank you for your time and help!

[Torsten]

As far as I can see from your code, solutions of previous time steps are never invoked.

RKCalcDyn3D works with the actual values for t and y.