33 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Angie vor etwa 19 Stunden
Bearbeitet: Torsten vor etwa 10 Stunden
In MATLAB Online öffnen
Hello everyone,
I am trying to solve some differential equation in matlab using ode45. The equations are as shown below:
I have written the equations in state space form where ydot_1(i) is , ydot_2(i) is
, ydot_3(i) is
and ydot_4(i) is
. In the equation of motion of x, eq.2, x terms are not in terms of “ i “ but I have written them like that (so ydot_3(i) is
and ydot_4(i) is
) because of the loop that I use on the function below. So the “for” loop, for each time step, solves first the ydot terms for each mode for i = 1:n and then gets out of the loop, goes to the new time step
and then solve again for each mode. However, since for each time step it will start from the first mode, and ydot_1(i) and ydot_2(i) depend on ydot_3(i) and ydot_4(i) , then I am solving the equation of x, eq.2, for each mode. For the new time step, before starting with the first mode again I am defining the
ydot_1=y(1:n); %displacement beam
ydot_2=y(n+1:2*n); %velocity beam
ydot_3=y(2*n+1:3*n); %displacement mass
ydot_4=y(3*n+1:4*n); %velocity mass
v_m=0;
x_m=0;
x_c_prev = x_c_prev_modes{1};
so it can start with the data from the previous first mode (first mode of the previous time step).
My problem is that when I solve for 1 mode only, so for n=1 and get the solution du_1= y(:,n+i).*sin(i*pi.* x_c_values_new); and then solve for n=2 but still get the solution du_1= y(:,n+i).*sin(i*pi.* x_c_values_new); for the first mode only, these first modes do not match. The solution have the same shape but differ in amplitude. The same happens with I use n=20, for example, but again plot only the first mode du_1. The solution changes from the first mode obtained when n=2 and from the first mode obtained for n=1. I have checked for days and cant see my mistake so have reached desperation.
Would really appreciate your help! Thank youuu!
clear all
clc
global xi xi_b r_w r_m mu_s mu_k v alpha P l a w l_0 n x_c_prev x_c_prev_modes
dt=0.01;
tspan = (0:dt:150).';
% parameters
xi=0.14;
xi_b=0.0025;
r_w=4.5;
r_m=4.8;
mu_s=0.5;
mu_k=0.2 ;
alpha=300;
P=0.002;
l=1/75;
l_0=0.25;
a=0;
w=9;
v=0.005;
% Initialize y_c_prev and the array to store x_c values
x_c_prev = 0;
global t_values x_c_values; % Declare t_values and x_c_values as global
t_values = tspan; % Store time values
x_c_values = zeros(size(tspan)); % Initialize array to store x_c values
n=1;
y0 = zeros(4*n,1);
y0(3*n+1) = -v;
% Initialize y_c_prev
x_c_prev = 0;
% Define your initial value for x_c_prev
x_c_prev_initial = 0;
% Initialize x_c_prev_modes, assuming you have 'n' modes
x_c_prev_modes = cell(1, n);
% Initialize the first mode's x_c_prev
x_c_prev_modes{1} = x_c_prev_initial;
[t,y] = ode45(@rhs_cont4,tspan,y0);
zero_indices = find(x_c_values == 0);
% Iterate over zero indices and replace with mean of neighboring non-zero elements
for i = 1:length(zero_indices)
idx = zero_indices(i);
left_idx = idx - 1;
right_idx = idx + 1;
% Find the nearest non-zero elements
while left_idx > 0 && x_c_values(left_idx) == 0
left_idx = left_idx - 1;
end
while right_idx <= length(x_c_values) && x_c_values(right_idx) == 0
right_idx = right_idx + 1;
end
% Replace zero with mean of neighboring non-zero elements
if left_idx > 0 && right_idx <= length(x_c_values)
x_c_values(idx) = (x_c_values(left_idx) + x_c_values(right_idx)) / 2;
elseif left_idx > 0
x_c_values(idx) = x_c_values(left_idx);
elseif right_idx <= length(x_c_values)
x_c_values(idx) = x_c_values(right_idx);
else
% If there are no non-zero neighbors, we can't replace it, so skip
continue;
end
end
x_c_values_new=x_c_values;
for i=1:1
u_1= y(:,i).*sin(i*pi.* x_c_values_new);
du_1= y(:,n+i).*sin(i*pi.* x_c_values_new);
end
figure
plot(t,du_1)
legend("beam vel mode 1")
and the function is saved separately
function ydot=rhs_cont4(t,y)
global xi xi_b r_w r_m mu_s mu_k v alpha P a w l_0 n x_c_prev x_c_prev_modes
ydot_1=y(1:n); %displacement beam
ydot_2=y(n+1:2*n); %velocity beam
ydot_3=y(2*n+1:3*n); %displacement mass
ydot_4=y(3*n+1:4*n); %velocity mass
v_m=0;
x_m=0;
x_c_prev = x_c_prev_modes{1};
for i = 1:n
x_m=y(i)*sin(i*pi*x_c_prev)+x_m;
x_c=y(2*n+i)+v*t+l_0-x_m;
x_c_prev = x_c; % Update x_c_prev for the next iteration
x_c_prev_modes{i} = x_c_prev;
%Store x_c values at the corresponding time step
global t_values x_c_values; % Access global variables
[~, idx] = min(abs(t - t_values)); % Find the index of the closest time step
x_c_values(idx) = x_c; % Store x_c at the corresponding time step
v_m=y(n+i)*sin(i*pi*x_c)+v_m;
ydot_1(i) = y(n+i);
ydot_2(i) = -2*xi_b*r_w*y(n+i)-(r_w*i)^2*y(i)+2*sign(y(3*n+i)+v-v_m)*(mu_k+(mu_s-mu_k)*exp(-alpha*abs(y(3*n+i)+v-v_m)))*P*sin(i*pi*x_c);
ydot_3(i) =y(3*n+i);
ydot_4(i) =-2*xi*y(3*n+i)-y(2*n+i)-r_m*sign(y(3*n+i)+v-v_m)*(mu_k+(mu_s-mu_k)*exp(-alpha*abs(y(3*n+i)+v-v_m)))*P+a*w^2*sin(w*t);
end
ydot = [ydot_1; ydot_2; ydot_3; ydot_4];
end
8 Kommentare 6 ältere Kommentare anzeigen6 ältere Kommentare ausblenden
6 ältere Kommentare anzeigen6 ältere Kommentare ausblenden
Angie vor etwa 19 Stunden
Direkter Link zu diesem Kommentar
https://de.mathworks.com/matlabcentral/answers/2135888-solving-coupled-differential-equations-using-ode45-for-n-modes#comment_3206973
to make it hopefully clearer, when i run for n=1 and n=2 i get the du_1 plots in the figure. so the one of n=1 is slighly bigger in amplitude (the difference increases when i include more modes). since i am plotting the first modes, i think i should get the same results.
I have printed before the results per timestep, for n=1 and i=1 and then for n=2 but for i=1 only. intialy they are the same but then change after some time. if i look at the code though, i dont see where i am making the mistake.
Torsten vor etwa 18 Stunden
Direkter Link zu diesem Kommentar
https://de.mathworks.com/matlabcentral/answers/2135888-solving-coupled-differential-equations-using-ode45-for-n-modes#comment_3207048
Is there a publication from which you took your equations ?
Angie vor etwa 17 Stunden
Direkter Link zu diesem Kommentar
https://de.mathworks.com/matlabcentral/answers/2135888-solving-coupled-differential-equations-using-ode45-for-n-modes#comment_3207078
- stick slip vibration of a moving osc on an axillary flex beam.pdf
Yes, I am attaching it here. Equation 1 in my post is eq.33 in the paper, eq.2 is eq.38, eq.3 - eq.5 are eq.34 - eq.37. Equation 6 in the post is eq.20 in the paper.
Torsten vor etwa 14 Stunden
Direkter Link zu diesem Kommentar
https://de.mathworks.com/matlabcentral/answers/2135888-solving-coupled-differential-equations-using-ode45-for-n-modes#comment_3207148
Why don't you follow the article and set up system (40) ?
Angie vor etwa 14 Stunden
Direkter Link zu diesem Kommentar
https://de.mathworks.com/matlabcentral/answers/2135888-solving-coupled-differential-equations-using-ode45-for-n-modes#comment_3207173
Bearbeitet: Angie vor etwa 14 Stunden
In the system (40), the equations I wrote are expressed in matrix form. There are N+1 second-order differential equations (and 1 for x_c): N equations for T''_i and one equation for x'' . When these are written in state-space form, we have 2N+2 equations. In my approach, I solve for 4N equations, where the first 2N pertain to T_i and the last 2N to x_i. I solve the equation for x for each mode because the equations are coupled. To find the solutions, I use a for loop to iterate over each N and employ `ode45`. I m not sure if there is another way to do it, but i am indeed solving the equaitons in (40).
Torsten vor etwa 14 Stunden
Direkter Link zu diesem Kommentar
https://de.mathworks.com/matlabcentral/answers/2135888-solving-coupled-differential-equations-using-ode45-for-n-modes#comment_3207183
Bearbeitet: Torsten vor etwa 14 Stunden
I m not sure if there is another way to do it, but i am indeed solving the equaitons in (40).
You don't. The number of equations to be solved should be 1 algebraic equation for x_c, 2 first-order ordinary differential equations for xi and xidot and 2*N first-order ordinary differential equations for T1,...,TN,T1dot,...TNdot. In total 2*N+3 equations from which one is algebraic and 2*N+2 are differential.
Angie vor etwa 13 Stunden
Direkter Link zu diesem Kommentar
https://de.mathworks.com/matlabcentral/answers/2135888-solving-coupled-differential-equations-using-ode45-for-n-modes#comment_3207203
Bearbeitet: Angie vor etwa 13 Stunden
In MATLAB Online öffnen
intially i solved like that so for 2*N+2 differential equation and had the same problem. The reason i changed to 4*N differential equation is because the equations of T_i depend on x. Starting at the first time step, when i enter the for loop and solve for i=1:n, i have the solutions of T per each mode but only one solution of x (where the f term in the equation of x depends of T_1 to T_N). When i go out of the loop and move to the new time step but start with i=1, i wanted the solution of x inside T to be of the first mode only (and not have the acumulation of all the modes from the first time step), since i was solving T for the first mode again (new time step). that is why i started to solve x for each mode as well.
when i solved for 2*N+2, i had almost the same function but with ydot_3 and ydot_4 not depending on i, as shown below.
function ydot=rhs_cont3(t,y)
global xi xi_b r_w r_m mu_s mu_k v alpha P l a w l_0 n x_c_prev x_c_prev_modes x_m v_m
ydot_1=y(1:n); %displacement beam
ydot_2=y(n+1:2*n); %velocity beam
ydot_3=y(2*n+1); %displacement mass
ydot_4=y(2*n+2); %velocity mass
v_m=0;
x_m=0;
x_c_prev = x_c_prev_modes{1};
for i = 1:n
x_m=y(i)*sin(i*pi*x_c_prev)+x_m;
x_c=y(2*n+1)+v*t+l_0-x_m;
x_c_prev = x_c; % Update y_c_prev for the next iteration
x_c_prev_modes{i} = x_c_prev;
%Store y_c values at the corresponding time step
global t_values x_c_values; % Access global variables
[~, idx] = min(abs(t - t_values)); % Find the index of the closest time step
x_c_values(idx) = x_c; % Store y_c at the corresponding time step
v_m=y(n+i)*sin(i*pi*x_c)+v_m;
ydot_1(i) = y(n+i);
ydot_2(i) = -2*xi_b*r_w*y(n+i)-(r_w*i)^2*y(i)+2*sign(y(2*n+2)+v-v_m)*(mu_k+(mu_s-mu_k)*exp(-alpha*abs(y(2*n+2)+v-v_m)))*P*sin(i*pi*x_c);
ydot_3 =y(2*n+2);
ydot_4 =-2*xi*y(2*n+2)-y(2*n+1)-r_m*sign(y(2*n+2)+v-v_m)*(mu_k+(mu_s-mu_k)*exp(-alpha*abs(y(2*n+2)+v-v_m)))*P+a*w^2*sin(w*t);
% if i == 1
% fprintf('Time: %.4f, ydot_4: %.7f\n ', t, ydot_4);
% end
end
ydot = [ydot_1; ydot_2; ydot_3; ydot_4];
and in the main function the initial conditions are defined for 2*N+2
n=1;
y0 = zeros(2*n+2,1);
y0(2*n+2) = -v;
the rest of the code is pretty much the same. in this case, the solutions seem further away from the ones in the paper.
Torsten vor etwa 11 Stunden
Direkter Link zu diesem Kommentar
https://de.mathworks.com/matlabcentral/answers/2135888-solving-coupled-differential-equations-using-ode45-for-n-modes#comment_3207248
Bearbeitet: Torsten vor etwa 10 Stunden
intially i solved like that so for 2*N+2 differential equation and had the same problem.
ode45 must solve for 2*N+3 equations where one of them is algebraic. Look at the mass matrix in equation (40). x_c is an algebraic solution variable of ode45, and the respective row of the mass matrix is set to 0. No need to use x_c_prev or find or min or whatever you use in your code to cope with and determine x_c.
When i go out of the loop and move to the new time step but start with i=1, i wanted the solution of x inside T to be of the first mode only (and not have the acumulation of all the modes from the first time step), since i was solving T for the first mode again (new time step). that is why i started to solve x for each mode as well.
I don't know what you mean here. Equation (35) cleary shows how the f_i are computed using x_c, xi, xidot, T1,...,TN,T1dot,...,TNdot. And no more than these 2*N+3 variables are needed.
Melden Sie sich an, um zu kommentieren.
Melden Sie sich an, um diese Frage zu beantworten.
Antworten (0)
Melden Sie sich an, um diese Frage zu beantworten.
Es ist ein Fehler aufgetreten
Da Änderungen an der Seite vorgenommen wurden, kann diese Aktion nicht abgeschlossen werden. Laden Sie die Seite neu, um sie im aktualisierten Zustand anzuzeigen.
Website auswählen
Wählen Sie eine Website aus, um übersetzte Inhalte (sofern verfügbar) sowie lokale Veranstaltungen und Angebote anzuzeigen. Auf der Grundlage Ihres Standorts empfehlen wir Ihnen die folgende Auswahl: .
Sie können auch eine Website aus der folgenden Liste auswählen:
Amerika
- América Latina (Español)
- Canada (English)
- United States (English)
Europa
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- Deutsch
- English
- Français
- United Kingdom(English)
Asien-Pazifik
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)
Kontakt zu Ihrer lokalen Niederlassung