HOW TO: Using each step result in subsequent step during integration via ODE

Hello, I am trying to solve these equation of motion using ode45. The implementation process in Matrix format is presented below.

x0=zeros(16,1);
opts=odeset('RelTol',1e-6);

% Freq = linspace(1,4000,50);

ti = 0;   % initial time step
tf = 1; % final time step
dt = 0.001 ;            % size of the time step
nt = fix((tf-ti)/dt)+1;   % number of time steps
tau = 0:dt:tf;

    W = 4000
    [t, Xx] = ode45(@gearCH, tau, x0, opts, W);
    

function Xx=gearCH(t,x0,W)

m1 = 2.8;
i1 = 0.003655;
j1 = 0.00731;

m2 = m1;
i2 = i1;
j2 = j1;

kby1 = 5e9;
kbz1 = 4e8;
kb0 = 3.23e6;

kby2 = kby1;
kbz2 = kbz1;

r1 = 70e-3;
r2 = r1;

beta=deg2rad(25.232);

e5= .01e-6;

torload = 450;      % Output Torque (lb-in, N-m)


et = e5 * sin(W * t);

M = blkdiag(m1, m1, i1/(r1^2), j1/(r1^2), m2, m2, i2/(r2^2), j2/(r2^2)); %mass
Kb = blkdiag (kby1, kbz1, kb0/(r1^2), 0, kby2, kbz2, kb0/(r2^2), 0); %bearing stiff
Q = [cos(beta) -sin(beta) sin(beta) cos(beta) -cos(beta) sin(beta) sin(beta) cos(beta)]'; %Q

Kt = 3.6e8; %%use as time invariant here

Km = Kt *( Q* Q');


Fb = [0 0 0 torload/r1 0 0 0 torload/r2]';
Fc = (Kt*et)*Q;



cma = (2*0.03) * sqrt(Kt /((1/r1) + (1/r2)));
cG = cma*( Q* Q');  %gear estimated damping
cB = 2.5e-5*Kb; %bearing estimated damping


K = Kb + Km; %stiffness matrix
C = cG + cB; %damping matrix

P = Fb + Fc; %force

Xx=[x0(9:16);inv(M)*(P-C*x0(9:16)-K*x0(1:8))];
end

As can be seen, this is done using ode45 and also, I have not use p(t) and p(t) dot, to multiply K and C as required. How can I get the displacements X, for each time so I can calulate p(t), multiply it by K and C for the next integration during the whole process? It seems I have to implement Rung-Kutta manually, maybe not. But how can I achieve that given this 8 second order odes? Thank you

Answer 1

[Joe]

I am not sure this is doable using matlab. The equations are not common, at least to me. Try SciPy maybe

Answer 2

[Torsten]

All variables needed to compute p(t) are part of the solution vector x which is input to "gearCH". So I don't see the problem to compute p(t) from them.

[Presley]

The variables seem confusing and I have modified them. I intend to get Xx result at time t and use it to calcalute p, then multiply it by K to get the new K. This p(t) is also differentiated and multiplied with C to get the new C. Then solve the equation again for the next integration (t+1). This is repeated until final step.

As it is, ode just solves and give all Xx for all time. I need Xx(t) value to get the next Xx(t+1) output

[Torsten]

In "gearCH", you get the X(t) value, and p(t) can be computed as

pt = (x0(1)-x0(5)+x0(4)+x0(8))*cos(beta) + (-x0(2)+x0(6)+x0(3)+x0(7))*sin(beta) - et

and its derivative dp/dt as

ptt = (x0(8+1)-x0(8+5)+x0(8+4)+x0(8+8))*cos(beta) + (-x0(8+2)+x0(8+6)+x0(8+3)+x0(8+7))*sin(beta) - e5*W*cos(W*t)

After you supplied Xx, ode45 will do the next time step.

[Presley]

Thanks for your response. However, when I checked back I realised I need not have to use pt and its derivative when solving in the matrix form. So, running the code as I had initially posted above using ode45 in matlab, and plotting time vs Xx(:,16 ) gives the figure below.

Now, I tried to use a manual implementation of runge-kutta 4th order to check the above since the intention was to use a higher order for this problem as required. But the plot is not looking good. The code and figure below:

    %% rungee kutta 4th implementation check

clear all; clc; close all;

x0=zeros(16,1);
h = 5e-2;
tEnd = 5;
coeff_1per6 = 1/6;
    W = 5000;
yHist = [];
tHist = [];

storeID = 0;
tic
for tStep = 0:h:tEnd
%
k1 = h*gearCH(tStep,x0,W);
k2 = h*gearCH(tStep+0.5*h,x0+0.5*h*k1,W);
k3 = h*gearCH(tStep+0.5*h,x0+0.5*h*k2,W);
k4 = h*gearCH(tStep+h,x0+h*k3,W);
y1 = x0+coeff_1per6*(k1+2*k2+2*k3+k4);


    storeID = storeID+1;
    yHist(1:16,storeID) = y1;
    tHist(storeID) = tStep;

x0 = y1;

end
toc

plot(tHist,yHist(16,:))
    

Any help is appreciated. Thanks

[Torsten]

RK4 is of higher order than ode45 ? 4 is bigger than 5 ?

Try a smaller time step with your self-written code. Maybe this can prevent it from diverging.

I cannot understand why you think a self-written numerical integration code could be helpful here, especially because you use the same function routine "gearCH" that supplies the same time derivatives for both codes.

[Presley]

The author of the text used runge-kutta-vener 6th order instead of rung-kutta fehlberg. Now, when I used ode45, the subsequent post processing of results (not the wave plot above) does not yield acuurate results as the text. So, I was manually trying RK4 to check if it works before proceeding to veners 6th order..see the >> Rk4 vern.. I will implement it directly now and give feedback. thanks

[Walter Roberson]

Did you happen to experiment with ode78 or ode89 ?

[Torsten]

You will get exact solutions also with a low order code if you strengthen the tolerances in the options setting. So I suggest using ODE45 with smaller values for RelTol and/or AbsTol.

Answer 3

[Sam Chak]

Hi @Presley

From your provided formula, both and can be directly computed from the ode45 solution vector x. Could you please verify if they have been plotted correctly?

% Freq = linspace(1, 4000, 50);
ti     = 0;                     % initial time step
tf     = 0.02;                  % final time step
dt     = 1/20000;               % size of the time step
% nt   = fix((tf-ti)/dt)+1;     % number of time steps

% call ode45 solver
tspan  = 0:dt:tf;
x0     = zeros(16, 1);
opts   = odeset('RelTol', 1e-6);
W      = 4000;
[t, x] = ode45(@gearCH, tspan, x0, opts, W);

% plot pt and dp/dt
beta   = deg2rad(25.232);
e5     = .01e-6;
et     = e5*sin(W*t);
pt     = (x(:,1) -  x(:,5) +  x(:,4) +  x(:,8))*cos(beta) + (-  x(:,2) +  x(:,6) +  x(:,3) +  x(:,7))*sin(beta) - et;
dpdt   = (x(:,9) - x(:,13) + x(:,12) + x(:,16))*cos(beta) + (- x(:,10) + x(:,14) + x(:,11) + x(:,15))*sin(beta) - e5*W*cos(W*t);
figure(1)
subplot(211)
plot(t, pt),   grid on, xlabel('t'), ylabel('p_{t}')
subplot(212)
plot(t, dpdt), grid on, xlabel('t'), ylabel('dp/dt')

% plot solution
figure(2)
for j = 1:16
    subplot(4,4,j)
    plot(t, x(:,j)), grid on
    title("x"+string(j));
end

% Equations of Motion
function dxdt = gearCH(t, x, W)

    % parameters
    m1      = 2.8;
    i1      = 0.003655;
    j1      = 0.00731;    
    m2      = m1;
    i2      = i1;
    j2      = j1;    
    kby1    = 5e9;
    kbz1    = 4e8;
    kb0     = 3.23e6;    
    kby2    = kby1;
    kbz2    = kbz1;    
    r1      = 70e-3;
    r2      = r1;    
    beta    = deg2rad(25.232);    
    e5      = .01e-6;    
    torload = 450;      % Output Torque (lb-in, N-m)
    et      = e5*sin(W*t);
    
    pt     = (x(1) -  x(5) +  x(4) +  x(8))*cos(beta) + (-  x(2) +  x(6) +  x(3) +  x(7))*sin(beta) - et;
    dpdt   = (x(9) - x(13) + x(12) + x(16))*cos(beta) + (- x(10) + x(14) + x(11) + x(15))*sin(beta) - e5*W*cos(W*t);
        
    % Mass matrix
    M       = blkdiag(m1, m1, i1/(r1^2), j1/(r1^2), m2, m2, i2/(r2^2), j2/(r2^2));
    % Bearing stiffness matrix
    Kb      = blkdiag (kby1, kbz1, kb0/(r1^2), 0, kby2, kbz2, kb0/(r2^2), 0);
    % Q array
    Q       = [cos(beta) -sin(beta)  sin(beta)  cos(beta) -cos(beta)  sin(beta)  sin(beta)  cos(beta)]';
        
    % unnamed parameter used in the computation of Km
    Kt      = 3.6e8;            % used as time invariant here
    % untitled Q-based square matrix
    Km      = pt*Kt*(Q*Q');     % <-- pt is injected here
    
    % untitled b-force
    Fb      = [0 0 0 torload/r1 0 0 0 torload/r2]';
    % untitled c-force
    Fc      = (Kt*et)*Q;
    
    % unnamed parameter used in the computation of cG
    cma     = 2*0.03*sqrt(Kt/(1/r1 + 1/r2));
    % gear-estimated damping marix
    cG      = dpdt*cma*(Q*Q');  % <-- dp/dt is injected here
    % bearing-estimated damping matrix
    cB      = (2.5e-5)*Kb;
    
    % True stiffness matrix
    K       = Kb + Km; 
    % True damping matrix
    C       = cG + cB; 
    
    % Total force
    F       = Fb + Fc;
    
    % Equations of Motion
    dxdt       = zeros(16, 1);
    dxdt(1:8)  = x(9:16);                       % kinematics
    dxdt(9:16) = M\(F - C*x(9:16) - K*x(1:8));  % dynamics
end

[Presley]

Thanks for the response, however, the problem continued under Torsten post. I was trying a manual implementation of the these solvers in a higher order. I will give feedback after trying the 6th order, Thanks