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

3 views (last 30 days)
Presley on 31 Jul 2023
Commented: Torsten on 7 Oct 2023
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;
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');
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

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

Torsten on 31 Jul 2023
Edited: Torsten on 31 Jul 2023
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.
Walter Roberson on 7 Oct 2023
Did you happen to experiment with ode78 or ode89 ?
Torsten on 7 Oct 2023
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.

Sam Chak on 7 Oct 2023
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
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;
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
% 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 on 7 Oct 2023
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

### Categories

Find more on Programming in Help Center and File Exchange

R2022b

### Community Treasure Hunt

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

Start Hunting!