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

Question

Presley 2023-7-31

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2002617-how-to-using-each-step-result-in-subsequent-step-during-integration-via-ode

评论： Torsten 2023-10-7

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Joe 2023-7-31

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2002617-how-to-using-each-step-result-in-subsequent-step-during-integration-via-ode#answer_1281027

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Torsten 2023-7-31

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2002617-how-to-using-each-step-result-in-subsequent-step-during-integration-via-ode#answer_1281097

编辑：Torsten 2023-7-31

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.

7 个评论
显示 5更早的评论隐藏 5更早的评论

Presley 2023-10-6

在 MATLAB Online 中打开

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

Walter Roberson 2023-10-7

Did you happen to experiment with ode78 or ode89 ?

Torsten 2023-10-7

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 2023-10-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2002617-how-to-using-each-step-result-in-subsequent-step-during-integration-via-ode#answer_1328344

在 MATLAB Online 中打开

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Presley 2023-10-7

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

请先登录，再进行评论。