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

3 views (last 30 days)

Show older comments

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

##### 0 Comments

### Answers (3)

Joe
on 31 Jul 2023

##### 0 Comments

Torsten
on 31 Jul 2023

Edited: Torsten
on 31 Jul 2023

##### 7 Comments

Torsten
on 7 Oct 2023

Sam Chak
on 7 Oct 2023

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

### See Also

### Community Treasure Hunt

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

Start Hunting!