solving coupled differential equations using ode45 for “n” modes
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];
endHello 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 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 ode45, differential equations, for loop MATLAB Answers — New Questions
