How to solve 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 )b ecause of the loop that I use on the function below. So the “for” loop, for each time step, lets say starting at 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 again, 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);
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")
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
ydot_1=zeros(n,1); %displacement beam
ydot_2=zeros(n,1); %velocity beam
ydot_3=zeros(n,1); %displacement mass
ydot_4=zeros(n,1); %velocity mass%
v_m=0;
x_m=0;
x_c_prev = x_c_prev_modes{1};
% fprintf(‘x_c_prev %.7fn ‘, x_c_prev);
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 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(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);
%
% if i == 1
% fprintf(‘Time: %.4f, ydot_1: %.7fn, ydot_2: %.7fn, ydot_3: %.7fn, ydot_4: %.7fn ‘, t, ydot_1(i), ydot_2(i), ydot_3(i), ydot_4(i));
% end
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 )b ecause of the loop that I use on the function below. So the “for” loop, for each time step, lets say starting at 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 again, 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);
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")
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
ydot_1=zeros(n,1); %displacement beam
ydot_2=zeros(n,1); %velocity beam
ydot_3=zeros(n,1); %displacement mass
ydot_4=zeros(n,1); %velocity mass%
v_m=0;
x_m=0;
x_c_prev = x_c_prev_modes{1};
% fprintf(‘x_c_prev %.7fn ‘, x_c_prev);
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 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(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);
%
% if i == 1
% fprintf(‘Time: %.4f, ydot_1: %.7fn, ydot_2: %.7fn, ydot_3: %.7fn, ydot_4: %.7fn ‘, t, ydot_1(i), ydot_2(i), ydot_3(i), ydot_4(i));
% end
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 )b ecause of the loop that I use on the function below. So the “for” loop, for each time step, lets say starting at 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 again, 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);
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")
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
ydot_1=zeros(n,1); %displacement beam
ydot_2=zeros(n,1); %velocity beam
ydot_3=zeros(n,1); %displacement mass
ydot_4=zeros(n,1); %velocity mass%
v_m=0;
x_m=0;
x_c_prev = x_c_prev_modes{1};
% fprintf(‘x_c_prev %.7fn ‘, x_c_prev);
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 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(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);
%
% if i == 1
% fprintf(‘Time: %.4f, ydot_1: %.7fn, ydot_2: %.7fn, ydot_3: %.7fn, ydot_4: %.7fn ‘, t, ydot_1(i), ydot_2(i), ydot_3(i), ydot_4(i));
% end
end
ydot = [ydot_1; ydot_2; ydot_3; ydot_4];
end ode45, differential equations MATLAB Answers — New Questions
