Tag Archives: matlab
maximum Lyapunov exponent diagram
hello, I’m working on discrete dynamical system in 2 dimesntion
t’m trying to plot a digram a paramete versus the maximum lyapunov exponent, i searched more about it but didn’t reach to anything.
any one have idea or could help me and thank you.hello, I’m working on discrete dynamical system in 2 dimesntion
t’m trying to plot a digram a paramete versus the maximum lyapunov exponent, i searched more about it but didn’t reach to anything.
any one have idea or could help me and thank you. hello, I’m working on discrete dynamical system in 2 dimesntion
t’m trying to plot a digram a paramete versus the maximum lyapunov exponent, i searched more about it but didn’t reach to anything.
any one have idea or could help me and thank you. dynamical system MATLAB Answers — New Questions
Timestep stability in a 1D heat diffusion model
Hi,
I have a 1D heat diffusion code which I was using on a timescale of 10s of years and I am now trying to use the same code to work on a scale of millions of years. Obviously if I keep my timestep the same this will take ages to calculate but if I increase my timestep I encounter numerical stability issues.
My questions are:
How should I approach this problem?
What affects the maximum stable timestep? And how do I calculate this?
Many thanks,
Alex
close all
clear all
dx = 4; % discretization step in m
dt = 0.0000001; % timestep in Myrs
h=1000; % height of box in m
nx=h/dx+1;
model_lenth=1; %length of model in Myrs
nt=ceil(model_lenth/dt)+1; % number of tsteps to reach end of model
kappa = 1e-6; % thermal diffusivity
x=0:dx:0+h; % finite difference mesh
T=38+0.05.*x; % initial T=Tm everywhere …
time=zeros(1,nt);
t=0;
Tnew = zeros(1,nx);
%Lower sill
sill_1_thickness=18;
Sill_1_top_position=590;
Sill_1_top=ceil(Sill_1_top_position/dx);
Sill_1_bottom=ceil((Sill_1_top_position+sill_1_thickness)/dx);
%Upper sill
sill_2_thickness=10;
Sill_2_top_position=260;
Sill_2_top=ceil(Sill_2_top_position/dx);
Sill_2_bottom=ceil((Sill_2_top_position+sill_2_thickness)/dx);
%Temperature of dolerite intrusions
Tm=1300;
T(Sill_1_top:Sill_1_bottom)=Tm; %Apply temperature to intrusion 1
% unit conversion to SI:
secinmyr=24*3600*365*1000000; % dt in sec
dt=dt*secinmyr;
%Plot initial conditions
figure(1), clf
f1 = figure(1); %Make full screen
set(f1,’Units’, ‘Normalized’, ‘OuterPosition’, [0 0 1 1]);
plot (T,x,’LineWidth’,2)
xlabel(‘T [^oC]’)
ylabel(‘x[m]’)
axis([0 1310 0 1000])
title(‘ Initial Conditions’)
set(gca,’YDir’,’reverse’);
%Main calculation
for it=1:nt
%Apply temperature to upper intrusion
if it==10;
T(Sill_2_top:Sill_2_bottom)=Tm;
end
for i = 2:nx-1
Tnew(i) = T(i) + kappa*dt*(T(i+1) – 2*T(i) + T(i-1))/dx/dx;
end
Tnew(1) = T(1);
Tnew(nx) = T(nx);
time(it) = t;
T = Tnew; %Set old Temp to = new temp for next loop
tmyears=(t/secinmyr);
%Plot a figure which updates in the loop of temperature against depth
figure(2), clf
plot (T,x,’LineWidth’,2)
xlabel(‘T [^oC]’)
ylabel(‘x[m]’)
title([‘ Temperature against Depth after ‘,num2str(tmyears),’ Myrs’])
axis([0 1300 0 1000])
set(gca,’YDir’,’reverse’);%Reverse y axis
%Make full screen
f2 = figure(2);
set(f2,’Units’, ‘Normalized’, ‘OuterPosition’, [0 0 1 1]);
drawnow
t=t+dt;
endHi,
I have a 1D heat diffusion code which I was using on a timescale of 10s of years and I am now trying to use the same code to work on a scale of millions of years. Obviously if I keep my timestep the same this will take ages to calculate but if I increase my timestep I encounter numerical stability issues.
My questions are:
How should I approach this problem?
What affects the maximum stable timestep? And how do I calculate this?
Many thanks,
Alex
close all
clear all
dx = 4; % discretization step in m
dt = 0.0000001; % timestep in Myrs
h=1000; % height of box in m
nx=h/dx+1;
model_lenth=1; %length of model in Myrs
nt=ceil(model_lenth/dt)+1; % number of tsteps to reach end of model
kappa = 1e-6; % thermal diffusivity
x=0:dx:0+h; % finite difference mesh
T=38+0.05.*x; % initial T=Tm everywhere …
time=zeros(1,nt);
t=0;
Tnew = zeros(1,nx);
%Lower sill
sill_1_thickness=18;
Sill_1_top_position=590;
Sill_1_top=ceil(Sill_1_top_position/dx);
Sill_1_bottom=ceil((Sill_1_top_position+sill_1_thickness)/dx);
%Upper sill
sill_2_thickness=10;
Sill_2_top_position=260;
Sill_2_top=ceil(Sill_2_top_position/dx);
Sill_2_bottom=ceil((Sill_2_top_position+sill_2_thickness)/dx);
%Temperature of dolerite intrusions
Tm=1300;
T(Sill_1_top:Sill_1_bottom)=Tm; %Apply temperature to intrusion 1
% unit conversion to SI:
secinmyr=24*3600*365*1000000; % dt in sec
dt=dt*secinmyr;
%Plot initial conditions
figure(1), clf
f1 = figure(1); %Make full screen
set(f1,’Units’, ‘Normalized’, ‘OuterPosition’, [0 0 1 1]);
plot (T,x,’LineWidth’,2)
xlabel(‘T [^oC]’)
ylabel(‘x[m]’)
axis([0 1310 0 1000])
title(‘ Initial Conditions’)
set(gca,’YDir’,’reverse’);
%Main calculation
for it=1:nt
%Apply temperature to upper intrusion
if it==10;
T(Sill_2_top:Sill_2_bottom)=Tm;
end
for i = 2:nx-1
Tnew(i) = T(i) + kappa*dt*(T(i+1) – 2*T(i) + T(i-1))/dx/dx;
end
Tnew(1) = T(1);
Tnew(nx) = T(nx);
time(it) = t;
T = Tnew; %Set old Temp to = new temp for next loop
tmyears=(t/secinmyr);
%Plot a figure which updates in the loop of temperature against depth
figure(2), clf
plot (T,x,’LineWidth’,2)
xlabel(‘T [^oC]’)
ylabel(‘x[m]’)
title([‘ Temperature against Depth after ‘,num2str(tmyears),’ Myrs’])
axis([0 1300 0 1000])
set(gca,’YDir’,’reverse’);%Reverse y axis
%Make full screen
f2 = figure(2);
set(f2,’Units’, ‘Normalized’, ‘OuterPosition’, [0 0 1 1]);
drawnow
t=t+dt;
end Hi,
I have a 1D heat diffusion code which I was using on a timescale of 10s of years and I am now trying to use the same code to work on a scale of millions of years. Obviously if I keep my timestep the same this will take ages to calculate but if I increase my timestep I encounter numerical stability issues.
My questions are:
How should I approach this problem?
What affects the maximum stable timestep? And how do I calculate this?
Many thanks,
Alex
close all
clear all
dx = 4; % discretization step in m
dt = 0.0000001; % timestep in Myrs
h=1000; % height of box in m
nx=h/dx+1;
model_lenth=1; %length of model in Myrs
nt=ceil(model_lenth/dt)+1; % number of tsteps to reach end of model
kappa = 1e-6; % thermal diffusivity
x=0:dx:0+h; % finite difference mesh
T=38+0.05.*x; % initial T=Tm everywhere …
time=zeros(1,nt);
t=0;
Tnew = zeros(1,nx);
%Lower sill
sill_1_thickness=18;
Sill_1_top_position=590;
Sill_1_top=ceil(Sill_1_top_position/dx);
Sill_1_bottom=ceil((Sill_1_top_position+sill_1_thickness)/dx);
%Upper sill
sill_2_thickness=10;
Sill_2_top_position=260;
Sill_2_top=ceil(Sill_2_top_position/dx);
Sill_2_bottom=ceil((Sill_2_top_position+sill_2_thickness)/dx);
%Temperature of dolerite intrusions
Tm=1300;
T(Sill_1_top:Sill_1_bottom)=Tm; %Apply temperature to intrusion 1
% unit conversion to SI:
secinmyr=24*3600*365*1000000; % dt in sec
dt=dt*secinmyr;
%Plot initial conditions
figure(1), clf
f1 = figure(1); %Make full screen
set(f1,’Units’, ‘Normalized’, ‘OuterPosition’, [0 0 1 1]);
plot (T,x,’LineWidth’,2)
xlabel(‘T [^oC]’)
ylabel(‘x[m]’)
axis([0 1310 0 1000])
title(‘ Initial Conditions’)
set(gca,’YDir’,’reverse’);
%Main calculation
for it=1:nt
%Apply temperature to upper intrusion
if it==10;
T(Sill_2_top:Sill_2_bottom)=Tm;
end
for i = 2:nx-1
Tnew(i) = T(i) + kappa*dt*(T(i+1) – 2*T(i) + T(i-1))/dx/dx;
end
Tnew(1) = T(1);
Tnew(nx) = T(nx);
time(it) = t;
T = Tnew; %Set old Temp to = new temp for next loop
tmyears=(t/secinmyr);
%Plot a figure which updates in the loop of temperature against depth
figure(2), clf
plot (T,x,’LineWidth’,2)
xlabel(‘T [^oC]’)
ylabel(‘x[m]’)
title([‘ Temperature against Depth after ‘,num2str(tmyears),’ Myrs’])
axis([0 1300 0 1000])
set(gca,’YDir’,’reverse’);%Reverse y axis
%Make full screen
f2 = figure(2);
set(f2,’Units’, ‘Normalized’, ‘OuterPosition’, [0 0 1 1]);
drawnow
t=t+dt;
end time step; heat diffusion MATLAB Answers — New Questions
How to set the theorem of floquet for Linear variational equation with periodic equation?
I’m trying to analyze the stability regions of the equation of motion of a system with periodic coefficients (whirlflutter). I have only written its polynomial equation and examined its roots, but my plot did not match the reference implementation correctly.
If I wanna write the Hill function and use the ode45 code, is it possible for someone to correct my code and my approach to finding the monodromy matrix and calculating the Floquet multipliers?
clc;
clear all;
% Define parameters
N = 2; % Number of blades
I_thetaoverI_b = 2; % Moment of inertia pitch axis over I_b
I_psioverI_b = 2; % Moment of inertia yaw axis over I_b
C_thetaoverI_b = 0.00; % Damping coefficient over I_b
C_psioverI_b = 0.00; % Damping coefficient over I_b
h = 0.3; % rotor mast height, wing tip spar to rotor hub [m]
hoverR = 0.34;
R = h / hoverR; % radius [m]
gamma = 4; % lock number
V_knots = 325; % the rotor forward velocity [knots]
% Convert velocity from [knots] to [meters per second]
% 1 knot = 0.51444 meters per second
V = V_knots * 0.51444;
% Calculate angular velocity in radians per second
omega_rad_s = V / R;
% Convert angular velocity from radians per second to RPM
% 1 radian per second = (60 / (2 * pi)) RPM
Omega = omega_rad_s * (60 / (2 * pi));
freq_1_over_Omega = 1 / Omega;
%the flap moment aerodynamic coefficients for large V
M_b = -(1/10)*V;
M_u = 1/6;
%the propeller aerodynamic coefficients
H_u = V/2;
% Frequency ranges
f_pitch= 5:3:140;
f_yaw= 5:3:140;
divergence_map = [];
Rdivergence_map = [];
unstable = [];
% Modify the loop to iterate over time points
for i = 1:length(f_pitch)
for j = 1:length(f_yaw)
phi_steps = linspace(0, pi, 100); % Time steps within one period
for phi = phi_steps
% Angular frequencies for the current time point
w_omega_pitch = 2 * pi* f_pitch(i);
w_omega_yaw = 2 * pi * f_yaw(j);
K_psi = (w_omega_pitch^2) * I_psioverI_b;
K_theta = (w_omega_yaw^2) * I_thetaoverI_b;
% Define inertia matrix [M]
M_matrix = [I_thetaoverI_b + 1 + cos(2*phi), -sin(2*phi);
-sin(2*phi), I_psioverI_b + 1 – cos(2*phi)];
% Define damping matrix [D]
D11 = C_thetaoverI_b + h^2*gamma*H_u*(1 – cos(2*phi)) – gamma*M_b*(1 + cos(2*phi)) – (2 + 2*h*gamma*M_u)*sin(2*phi);
D12 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) – 2*(1 + cos(2*phi)) – 2*h*gamma*M_u*cos(2*phi);
D21 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) + 2*(1 – cos(2*phi)) – 2*h*gamma*M_u*cos(2*phi);
D22 = C_psioverI_b + h^2*gamma*H_u*(1 + cos(2*phi)) – gamma*M_b*(1 – cos(2*phi)) + (2 + 2*h*gamma*M_u)*sin(2*phi);
D_matrix = [D11, D12;
D21, D22];
% Define stiffness matrix [K]
K11 = K_theta – h*gamma*V*H_u*(1 – cos(2*phi)) + gamma*V*M_u*sin(2*phi);
K12 = -h*V*gamma*H_u*sin(2*phi) + gamma*V*M_u*(1 + cos(2*phi));
K21 = -h*gamma*V*H_u*sin(2*phi) – gamma*V*M_u*(1 – cos(2*phi));
K22 = K_psi – h*gamma*V*H_u*(1 + cos(2*phi)) – gamma*V*M_u*sin(2*phi);
K_matrix = [K11, K12;
K21, K22];
% Compute the system matrices
M11 = M_matrix(1, 1); M12 = M_matrix(1, 2); M21 = M_matrix(2, 1); M22 = M_matrix(2, 2);
D11 = D_matrix(1, 1); D12 = D_matrix(1, 2); D21 = D_matrix(2, 1); D22 = D_matrix(2, 2);
K11 = K_matrix(1, 1); K12 = K_matrix(1, 2); K21 = K_matrix(2, 1); K22 = K_matrix(2, 2);
% Find the roots of the polynomial equation
P0 = M11*M22-M12*M21;
P1 = (- D11*M22*1j – D22*M11*1j + M12*D21*j + D12*M21*j);
P2 = (D11*D22*(1j)^2 – K22*M11 – K11*M22 – D12*D21*(1j)^2 + M12*K21 + M21*K12);
P3 = (D11*K22*1j – D12*K21*1j – D21*K12*1j + D22*K11*1j);
P4 = K11*K22 – K12*K21;
P = roots([P0, P1, P2, P3, P4]);
r = 1 * P;
%Flutter
for k = 1:length(r)
if (real(r(k)) > 0)
if (imag(r(k)) <= 0)
unstable = [unstable; phi, K_psi, K_theta];
% Proximity check for 1/Ω divergence
if abs(real(r(k)) – freq_1_over_Omega) < 0.1
Rdivergence_map = [Rdivergence_map; phi, K_psi, K_theta];
end
end
end
end
%Divergence
if (real(det(K_matrix)) < 0)
divergence_map = [divergence_map; phi, K_psi, K_theta];
end
end
end
end
% Plotting section
figure;
hold on;
scatter(unstable(:,2), unstable(:,3), ‘filled’);
scatter(divergence_map(:,2), divergence_map(:,3), ‘filled’, ‘r’);
scatter(Rdivergence_map(:,2), Rdivergence_map(:,3), ‘filled’, ‘g’);
xlabel(‘K_psi’);
ylabel(‘K_theta’);
title(‘Whirl Flutter Diagram’);
legend(‘Flutter area’, ‘Divergence area’, ‘1/Ω Divergence area’);
hold off;I’m trying to analyze the stability regions of the equation of motion of a system with periodic coefficients (whirlflutter). I have only written its polynomial equation and examined its roots, but my plot did not match the reference implementation correctly.
If I wanna write the Hill function and use the ode45 code, is it possible for someone to correct my code and my approach to finding the monodromy matrix and calculating the Floquet multipliers?
clc;
clear all;
% Define parameters
N = 2; % Number of blades
I_thetaoverI_b = 2; % Moment of inertia pitch axis over I_b
I_psioverI_b = 2; % Moment of inertia yaw axis over I_b
C_thetaoverI_b = 0.00; % Damping coefficient over I_b
C_psioverI_b = 0.00; % Damping coefficient over I_b
h = 0.3; % rotor mast height, wing tip spar to rotor hub [m]
hoverR = 0.34;
R = h / hoverR; % radius [m]
gamma = 4; % lock number
V_knots = 325; % the rotor forward velocity [knots]
% Convert velocity from [knots] to [meters per second]
% 1 knot = 0.51444 meters per second
V = V_knots * 0.51444;
% Calculate angular velocity in radians per second
omega_rad_s = V / R;
% Convert angular velocity from radians per second to RPM
% 1 radian per second = (60 / (2 * pi)) RPM
Omega = omega_rad_s * (60 / (2 * pi));
freq_1_over_Omega = 1 / Omega;
%the flap moment aerodynamic coefficients for large V
M_b = -(1/10)*V;
M_u = 1/6;
%the propeller aerodynamic coefficients
H_u = V/2;
% Frequency ranges
f_pitch= 5:3:140;
f_yaw= 5:3:140;
divergence_map = [];
Rdivergence_map = [];
unstable = [];
% Modify the loop to iterate over time points
for i = 1:length(f_pitch)
for j = 1:length(f_yaw)
phi_steps = linspace(0, pi, 100); % Time steps within one period
for phi = phi_steps
% Angular frequencies for the current time point
w_omega_pitch = 2 * pi* f_pitch(i);
w_omega_yaw = 2 * pi * f_yaw(j);
K_psi = (w_omega_pitch^2) * I_psioverI_b;
K_theta = (w_omega_yaw^2) * I_thetaoverI_b;
% Define inertia matrix [M]
M_matrix = [I_thetaoverI_b + 1 + cos(2*phi), -sin(2*phi);
-sin(2*phi), I_psioverI_b + 1 – cos(2*phi)];
% Define damping matrix [D]
D11 = C_thetaoverI_b + h^2*gamma*H_u*(1 – cos(2*phi)) – gamma*M_b*(1 + cos(2*phi)) – (2 + 2*h*gamma*M_u)*sin(2*phi);
D12 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) – 2*(1 + cos(2*phi)) – 2*h*gamma*M_u*cos(2*phi);
D21 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) + 2*(1 – cos(2*phi)) – 2*h*gamma*M_u*cos(2*phi);
D22 = C_psioverI_b + h^2*gamma*H_u*(1 + cos(2*phi)) – gamma*M_b*(1 – cos(2*phi)) + (2 + 2*h*gamma*M_u)*sin(2*phi);
D_matrix = [D11, D12;
D21, D22];
% Define stiffness matrix [K]
K11 = K_theta – h*gamma*V*H_u*(1 – cos(2*phi)) + gamma*V*M_u*sin(2*phi);
K12 = -h*V*gamma*H_u*sin(2*phi) + gamma*V*M_u*(1 + cos(2*phi));
K21 = -h*gamma*V*H_u*sin(2*phi) – gamma*V*M_u*(1 – cos(2*phi));
K22 = K_psi – h*gamma*V*H_u*(1 + cos(2*phi)) – gamma*V*M_u*sin(2*phi);
K_matrix = [K11, K12;
K21, K22];
% Compute the system matrices
M11 = M_matrix(1, 1); M12 = M_matrix(1, 2); M21 = M_matrix(2, 1); M22 = M_matrix(2, 2);
D11 = D_matrix(1, 1); D12 = D_matrix(1, 2); D21 = D_matrix(2, 1); D22 = D_matrix(2, 2);
K11 = K_matrix(1, 1); K12 = K_matrix(1, 2); K21 = K_matrix(2, 1); K22 = K_matrix(2, 2);
% Find the roots of the polynomial equation
P0 = M11*M22-M12*M21;
P1 = (- D11*M22*1j – D22*M11*1j + M12*D21*j + D12*M21*j);
P2 = (D11*D22*(1j)^2 – K22*M11 – K11*M22 – D12*D21*(1j)^2 + M12*K21 + M21*K12);
P3 = (D11*K22*1j – D12*K21*1j – D21*K12*1j + D22*K11*1j);
P4 = K11*K22 – K12*K21;
P = roots([P0, P1, P2, P3, P4]);
r = 1 * P;
%Flutter
for k = 1:length(r)
if (real(r(k)) > 0)
if (imag(r(k)) <= 0)
unstable = [unstable; phi, K_psi, K_theta];
% Proximity check for 1/Ω divergence
if abs(real(r(k)) – freq_1_over_Omega) < 0.1
Rdivergence_map = [Rdivergence_map; phi, K_psi, K_theta];
end
end
end
end
%Divergence
if (real(det(K_matrix)) < 0)
divergence_map = [divergence_map; phi, K_psi, K_theta];
end
end
end
end
% Plotting section
figure;
hold on;
scatter(unstable(:,2), unstable(:,3), ‘filled’);
scatter(divergence_map(:,2), divergence_map(:,3), ‘filled’, ‘r’);
scatter(Rdivergence_map(:,2), Rdivergence_map(:,3), ‘filled’, ‘g’);
xlabel(‘K_psi’);
ylabel(‘K_theta’);
title(‘Whirl Flutter Diagram’);
legend(‘Flutter area’, ‘Divergence area’, ‘1/Ω Divergence area’);
hold off; I’m trying to analyze the stability regions of the equation of motion of a system with periodic coefficients (whirlflutter). I have only written its polynomial equation and examined its roots, but my plot did not match the reference implementation correctly.
If I wanna write the Hill function and use the ode45 code, is it possible for someone to correct my code and my approach to finding the monodromy matrix and calculating the Floquet multipliers?
clc;
clear all;
% Define parameters
N = 2; % Number of blades
I_thetaoverI_b = 2; % Moment of inertia pitch axis over I_b
I_psioverI_b = 2; % Moment of inertia yaw axis over I_b
C_thetaoverI_b = 0.00; % Damping coefficient over I_b
C_psioverI_b = 0.00; % Damping coefficient over I_b
h = 0.3; % rotor mast height, wing tip spar to rotor hub [m]
hoverR = 0.34;
R = h / hoverR; % radius [m]
gamma = 4; % lock number
V_knots = 325; % the rotor forward velocity [knots]
% Convert velocity from [knots] to [meters per second]
% 1 knot = 0.51444 meters per second
V = V_knots * 0.51444;
% Calculate angular velocity in radians per second
omega_rad_s = V / R;
% Convert angular velocity from radians per second to RPM
% 1 radian per second = (60 / (2 * pi)) RPM
Omega = omega_rad_s * (60 / (2 * pi));
freq_1_over_Omega = 1 / Omega;
%the flap moment aerodynamic coefficients for large V
M_b = -(1/10)*V;
M_u = 1/6;
%the propeller aerodynamic coefficients
H_u = V/2;
% Frequency ranges
f_pitch= 5:3:140;
f_yaw= 5:3:140;
divergence_map = [];
Rdivergence_map = [];
unstable = [];
% Modify the loop to iterate over time points
for i = 1:length(f_pitch)
for j = 1:length(f_yaw)
phi_steps = linspace(0, pi, 100); % Time steps within one period
for phi = phi_steps
% Angular frequencies for the current time point
w_omega_pitch = 2 * pi* f_pitch(i);
w_omega_yaw = 2 * pi * f_yaw(j);
K_psi = (w_omega_pitch^2) * I_psioverI_b;
K_theta = (w_omega_yaw^2) * I_thetaoverI_b;
% Define inertia matrix [M]
M_matrix = [I_thetaoverI_b + 1 + cos(2*phi), -sin(2*phi);
-sin(2*phi), I_psioverI_b + 1 – cos(2*phi)];
% Define damping matrix [D]
D11 = C_thetaoverI_b + h^2*gamma*H_u*(1 – cos(2*phi)) – gamma*M_b*(1 + cos(2*phi)) – (2 + 2*h*gamma*M_u)*sin(2*phi);
D12 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) – 2*(1 + cos(2*phi)) – 2*h*gamma*M_u*cos(2*phi);
D21 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) + 2*(1 – cos(2*phi)) – 2*h*gamma*M_u*cos(2*phi);
D22 = C_psioverI_b + h^2*gamma*H_u*(1 + cos(2*phi)) – gamma*M_b*(1 – cos(2*phi)) + (2 + 2*h*gamma*M_u)*sin(2*phi);
D_matrix = [D11, D12;
D21, D22];
% Define stiffness matrix [K]
K11 = K_theta – h*gamma*V*H_u*(1 – cos(2*phi)) + gamma*V*M_u*sin(2*phi);
K12 = -h*V*gamma*H_u*sin(2*phi) + gamma*V*M_u*(1 + cos(2*phi));
K21 = -h*gamma*V*H_u*sin(2*phi) – gamma*V*M_u*(1 – cos(2*phi));
K22 = K_psi – h*gamma*V*H_u*(1 + cos(2*phi)) – gamma*V*M_u*sin(2*phi);
K_matrix = [K11, K12;
K21, K22];
% Compute the system matrices
M11 = M_matrix(1, 1); M12 = M_matrix(1, 2); M21 = M_matrix(2, 1); M22 = M_matrix(2, 2);
D11 = D_matrix(1, 1); D12 = D_matrix(1, 2); D21 = D_matrix(2, 1); D22 = D_matrix(2, 2);
K11 = K_matrix(1, 1); K12 = K_matrix(1, 2); K21 = K_matrix(2, 1); K22 = K_matrix(2, 2);
% Find the roots of the polynomial equation
P0 = M11*M22-M12*M21;
P1 = (- D11*M22*1j – D22*M11*1j + M12*D21*j + D12*M21*j);
P2 = (D11*D22*(1j)^2 – K22*M11 – K11*M22 – D12*D21*(1j)^2 + M12*K21 + M21*K12);
P3 = (D11*K22*1j – D12*K21*1j – D21*K12*1j + D22*K11*1j);
P4 = K11*K22 – K12*K21;
P = roots([P0, P1, P2, P3, P4]);
r = 1 * P;
%Flutter
for k = 1:length(r)
if (real(r(k)) > 0)
if (imag(r(k)) <= 0)
unstable = [unstable; phi, K_psi, K_theta];
% Proximity check for 1/Ω divergence
if abs(real(r(k)) – freq_1_over_Omega) < 0.1
Rdivergence_map = [Rdivergence_map; phi, K_psi, K_theta];
end
end
end
end
%Divergence
if (real(det(K_matrix)) < 0)
divergence_map = [divergence_map; phi, K_psi, K_theta];
end
end
end
end
% Plotting section
figure;
hold on;
scatter(unstable(:,2), unstable(:,3), ‘filled’);
scatter(divergence_map(:,2), divergence_map(:,3), ‘filled’, ‘r’);
scatter(Rdivergence_map(:,2), Rdivergence_map(:,3), ‘filled’, ‘g’);
xlabel(‘K_psi’);
ylabel(‘K_theta’);
title(‘Whirl Flutter Diagram’);
legend(‘Flutter area’, ‘Divergence area’, ‘1/Ω Divergence area’);
hold off; stability theory, floquet theorem, characteristic multiplier, periodic coeffcients MATLAB Answers — New Questions
Reps in Friedman test
Hi, I want to do Friedman test on a test data where data from an individual has been collected for 18 days. I have 4 groups with 5, 7, 8, and 6 individuals in each group. If I am understanding correctly, then the data should be in the following format:
[D11 D11
D21 D21
D31 D31
D12 D12];
So basically , each column would have 18*n group individuals and each column would represent different group. I want to make sure I am understanding the reqs correctly. Also, how would I represent the reps in this case? The code gives me bad rep number.
Thank you.Hi, I want to do Friedman test on a test data where data from an individual has been collected for 18 days. I have 4 groups with 5, 7, 8, and 6 individuals in each group. If I am understanding correctly, then the data should be in the following format:
[D11 D11
D21 D21
D31 D31
D12 D12];
So basically , each column would have 18*n group individuals and each column would represent different group. I want to make sure I am understanding the reqs correctly. Also, how would I represent the reps in this case? The code gives me bad rep number.
Thank you. Hi, I want to do Friedman test on a test data where data from an individual has been collected for 18 days. I have 4 groups with 5, 7, 8, and 6 individuals in each group. If I am understanding correctly, then the data should be in the following format:
[D11 D11
D21 D21
D31 D31
D12 D12];
So basically , each column would have 18*n group individuals and each column would represent different group. I want to make sure I am understanding the reqs correctly. Also, how would I represent the reps in this case? The code gives me bad rep number.
Thank you. friedman test, non-parametric, statistics MATLAB Answers — New Questions
Combined Discrete-Continuous Simulation In SIMULINK
I am not even sure if that is the right title for this question.
Here’s the question:
How to implement the following dynamic system in Simulink?
Here’s the problem:
Time points are used to establish a uniform time grid on the interval , where the constant step size is . The solver is assumed to be at the beginning of the time step , where the state at the current time , is known, and the objective is to compute the next state through the relation , where is described as follows:
is first transformed to through some algebraic functions,
an IVP is set up with its initial condition set to . This IVP is then solved using a variable step ODE solver from to ,
the solution at time , is then mapped to the next state by some other algebraic mappings.
The following figure shows how the state is propagated thorugh time, in a single step:
In the output, I need the intermediate solutions (denoted by x) in each step, as well as the state variables .This can be easily done in a script, but I am looking for a way to do it in Simulink; which may be trivial problem, but I am really confused at this moment and I can’t really find a way to do it.
In some ways, in each step, the continuous-time subsystem waits for the input to be available, and once it is there, the external solver "waits" for the continuous-time subsystem to locally solve the continuous-time dynamical equations using an adaptive step-size solver, then the external solver receives this subsystem’s output and continuous to perform the remaining functions to produce . I feel like this should be logical workflow of the simulation, but I am not sure how Simulink treats problems of this type.I am not even sure if that is the right title for this question.
Here’s the question:
How to implement the following dynamic system in Simulink?
Here’s the problem:
Time points are used to establish a uniform time grid on the interval , where the constant step size is . The solver is assumed to be at the beginning of the time step , where the state at the current time , is known, and the objective is to compute the next state through the relation , where is described as follows:
is first transformed to through some algebraic functions,
an IVP is set up with its initial condition set to . This IVP is then solved using a variable step ODE solver from to ,
the solution at time , is then mapped to the next state by some other algebraic mappings.
The following figure shows how the state is propagated thorugh time, in a single step:
In the output, I need the intermediate solutions (denoted by x) in each step, as well as the state variables .This can be easily done in a script, but I am looking for a way to do it in Simulink; which may be trivial problem, but I am really confused at this moment and I can’t really find a way to do it.
In some ways, in each step, the continuous-time subsystem waits for the input to be available, and once it is there, the external solver "waits" for the continuous-time subsystem to locally solve the continuous-time dynamical equations using an adaptive step-size solver, then the external solver receives this subsystem’s output and continuous to perform the remaining functions to produce . I feel like this should be logical workflow of the simulation, but I am not sure how Simulink treats problems of this type. I am not even sure if that is the right title for this question.
Here’s the question:
How to implement the following dynamic system in Simulink?
Here’s the problem:
Time points are used to establish a uniform time grid on the interval , where the constant step size is . The solver is assumed to be at the beginning of the time step , where the state at the current time , is known, and the objective is to compute the next state through the relation , where is described as follows:
is first transformed to through some algebraic functions,
an IVP is set up with its initial condition set to . This IVP is then solved using a variable step ODE solver from to ,
the solution at time , is then mapped to the next state by some other algebraic mappings.
The following figure shows how the state is propagated thorugh time, in a single step:
In the output, I need the intermediate solutions (denoted by x) in each step, as well as the state variables .This can be easily done in a script, but I am looking for a way to do it in Simulink; which may be trivial problem, but I am really confused at this moment and I can’t really find a way to do it.
In some ways, in each step, the continuous-time subsystem waits for the input to be available, and once it is there, the external solver "waits" for the continuous-time subsystem to locally solve the continuous-time dynamical equations using an adaptive step-size solver, then the external solver receives this subsystem’s output and continuous to perform the remaining functions to produce . I feel like this should be logical workflow of the simulation, but I am not sure how Simulink treats problems of this type. simulink, simulation MATLAB Answers — New Questions
how can i make a parent – children hierarchy
Hello
i have two column one is number of hiarachy the second is the name of that hiarachy like this :
can you help make a code that will generate what is manually wrote in the third column , note that the file is large and contain long WBS path like this (2.2.3.5.4.5.7)
RegardsHello
i have two column one is number of hiarachy the second is the name of that hiarachy like this :
can you help make a code that will generate what is manually wrote in the third column , note that the file is large and contain long WBS path like this (2.2.3.5.4.5.7)
Regards Hello
i have two column one is number of hiarachy the second is the name of that hiarachy like this :
can you help make a code that will generate what is manually wrote in the third column , note that the file is large and contain long WBS path like this (2.2.3.5.4.5.7)
Regards data, excel, matlab, matrix, importing excel data, code MATLAB Answers — New Questions
How to calculate the voltage unbalance of a three phase signal?
I need to take a signal and do the symmetrical decomposition to find the negative and positive sequence.
I found Simulink’s Sequence Analyzer. It makes this decomposition of a signal.
But I would like a routine or function that would do this with my signal without me having to use simulink.I need to take a signal and do the symmetrical decomposition to find the negative and positive sequence.
I found Simulink’s Sequence Analyzer. It makes this decomposition of a signal.
But I would like a routine or function that would do this with my signal without me having to use simulink. I need to take a signal and do the symmetrical decomposition to find the negative and positive sequence.
I found Simulink’s Sequence Analyzer. It makes this decomposition of a signal.
But I would like a routine or function that would do this with my signal without me having to use simulink. negative sequence, positive sequence, function, matlab function, symmetrical decomposition MATLAB Answers — New Questions
Ho do I save the values of While at each hour?
For a while loop how I can calculate the values at each hour before do a break? As we know while loop only performs the condition statment and breaks the loop, then gives you the value at that specific hour only. I provide te example below :
for i=1:Time
while water_tank_soc(i-1) < Water_tank_capacity_max
m_net_o(i) =………
m_net_i(i) = ……
% Update SOC
water_tank_soc(i) = water_tank_soc(i-1) + m_net_i(i) – m_net_o(i);
if water_tank_soc(i) >= Water_tank_capacity_max
break
end
% Update for next time step
water_tank_soc(i-1) = water_tank_soc(i); % Update previous state for next iteration
end
endFor a while loop how I can calculate the values at each hour before do a break? As we know while loop only performs the condition statment and breaks the loop, then gives you the value at that specific hour only. I provide te example below :
for i=1:Time
while water_tank_soc(i-1) < Water_tank_capacity_max
m_net_o(i) =………
m_net_i(i) = ……
% Update SOC
water_tank_soc(i) = water_tank_soc(i-1) + m_net_i(i) – m_net_o(i);
if water_tank_soc(i) >= Water_tank_capacity_max
break
end
% Update for next time step
water_tank_soc(i-1) = water_tank_soc(i); % Update previous state for next iteration
end
end For a while loop how I can calculate the values at each hour before do a break? As we know while loop only performs the condition statment and breaks the loop, then gives you the value at that specific hour only. I provide te example below :
for i=1:Time
while water_tank_soc(i-1) < Water_tank_capacity_max
m_net_o(i) =………
m_net_i(i) = ……
% Update SOC
water_tank_soc(i) = water_tank_soc(i-1) + m_net_i(i) – m_net_o(i);
if water_tank_soc(i) >= Water_tank_capacity_max
break
end
% Update for next time step
water_tank_soc(i-1) = water_tank_soc(i); % Update previous state for next iteration
end
end save a values, while loop MATLAB Answers — New Questions
How to approach this type of question in matlab?
A ship travels on a straight line course described by y = (2005x)/6, where distances are measured in kilometers. The ship starts when x = -20 and ends when x = 40. Calculate the distance at closest approach to a lighthouse located at the coordinate origin (0,0). Do not solve this using a plot.A ship travels on a straight line course described by y = (2005x)/6, where distances are measured in kilometers. The ship starts when x = -20 and ends when x = 40. Calculate the distance at closest approach to a lighthouse located at the coordinate origin (0,0). Do not solve this using a plot. A ship travels on a straight line course described by y = (2005x)/6, where distances are measured in kilometers. The ship starts when x = -20 and ends when x = 40. Calculate the distance at closest approach to a lighthouse located at the coordinate origin (0,0). Do not solve this using a plot. matrix manipulation, homework, physics MATLAB Answers — New Questions
Understanding Gaussian Process Regression in Regression Learner App
Hi everyone,
I’m having some trouble understanding Gaussian Process Regression (GPR) options in the Regression Learner App. There are three main choices for GPR models:
Predefined Kernel: I can directly choose a kernel (Rational Quadratic, Squared Exponential, Matern 5/2, or Exponential) if I know which one suits my data best.
All GPR Models (non-optimizable): If I’m unsure which kernel to use, I can select this option to try all non-optimizable GPR models.
Optimizable GPR: This option allows the hyperparameters to be optimized and has even more Kernels available.
Regardless of whether I choose the optimizable or non-optimizable version, each kernel has hyperparameters that can be tuned.
Here are my questions:
Automatic Hyperparameter Selection: For the standard model without optimization, the kernel parameters (hyperparameters) are automatically selected and are an initial best guess rather than the optimal ones to my understanding. If so, how are they estimated/selected, and why aren’t they optimal?
Optimization Process: When the optimization option is selected, the hyperparameters are found by maximizing the Log Marginal Likelihood function. So, they are not guessed, but these are the hyperparameters that best describe the data. So finding the best value is basically the optimization process?
Choosing Non-Optimal Hyperparameters: Why would someone use the version without the best hyperparameters? One reason I’ve encountered is that optimization can take too long with the same data. Are there other reasons to choose non-optimizable models?
Best Practices: Are there standard practices or guidelines for when to use each version of the GPR models? Or one usually starts with non-optimizable version and then go to the optimizable one if the results are not sufficient enough?
Output Differences: Both versions output forecasted data and the covariance matrix. Does the optimized version provide any additional information or benefits?
Plotting Covariance: How can I plot the covariance that quantifies the uncertainty of my predictions? This is one of the main advantages of using GPR, and I want to visualize it. I have found some sample code here https://ch.mathworks.com/help/stats/gaussian-process-regression-models.html But I am a bit confused, because my input is a table actually, so I do not have only one predictor but a table. So far I have always used “predictedData = trainedModel.predictFcn(Table); but this gives me only the forecasted data without the 95% prediction intervals
I’m generally confused about the differences between both versions and the results they produce. Any insights or explanations would be greatly appreciated.
Thanks in advance!
Best regards,
GeorgiHi everyone,
I’m having some trouble understanding Gaussian Process Regression (GPR) options in the Regression Learner App. There are three main choices for GPR models:
Predefined Kernel: I can directly choose a kernel (Rational Quadratic, Squared Exponential, Matern 5/2, or Exponential) if I know which one suits my data best.
All GPR Models (non-optimizable): If I’m unsure which kernel to use, I can select this option to try all non-optimizable GPR models.
Optimizable GPR: This option allows the hyperparameters to be optimized and has even more Kernels available.
Regardless of whether I choose the optimizable or non-optimizable version, each kernel has hyperparameters that can be tuned.
Here are my questions:
Automatic Hyperparameter Selection: For the standard model without optimization, the kernel parameters (hyperparameters) are automatically selected and are an initial best guess rather than the optimal ones to my understanding. If so, how are they estimated/selected, and why aren’t they optimal?
Optimization Process: When the optimization option is selected, the hyperparameters are found by maximizing the Log Marginal Likelihood function. So, they are not guessed, but these are the hyperparameters that best describe the data. So finding the best value is basically the optimization process?
Choosing Non-Optimal Hyperparameters: Why would someone use the version without the best hyperparameters? One reason I’ve encountered is that optimization can take too long with the same data. Are there other reasons to choose non-optimizable models?
Best Practices: Are there standard practices or guidelines for when to use each version of the GPR models? Or one usually starts with non-optimizable version and then go to the optimizable one if the results are not sufficient enough?
Output Differences: Both versions output forecasted data and the covariance matrix. Does the optimized version provide any additional information or benefits?
Plotting Covariance: How can I plot the covariance that quantifies the uncertainty of my predictions? This is one of the main advantages of using GPR, and I want to visualize it. I have found some sample code here https://ch.mathworks.com/help/stats/gaussian-process-regression-models.html But I am a bit confused, because my input is a table actually, so I do not have only one predictor but a table. So far I have always used “predictedData = trainedModel.predictFcn(Table); but this gives me only the forecasted data without the 95% prediction intervals
I’m generally confused about the differences between both versions and the results they produce. Any insights or explanations would be greatly appreciated.
Thanks in advance!
Best regards,
Georgi Hi everyone,
I’m having some trouble understanding Gaussian Process Regression (GPR) options in the Regression Learner App. There are three main choices for GPR models:
Predefined Kernel: I can directly choose a kernel (Rational Quadratic, Squared Exponential, Matern 5/2, or Exponential) if I know which one suits my data best.
All GPR Models (non-optimizable): If I’m unsure which kernel to use, I can select this option to try all non-optimizable GPR models.
Optimizable GPR: This option allows the hyperparameters to be optimized and has even more Kernels available.
Regardless of whether I choose the optimizable or non-optimizable version, each kernel has hyperparameters that can be tuned.
Here are my questions:
Automatic Hyperparameter Selection: For the standard model without optimization, the kernel parameters (hyperparameters) are automatically selected and are an initial best guess rather than the optimal ones to my understanding. If so, how are they estimated/selected, and why aren’t they optimal?
Optimization Process: When the optimization option is selected, the hyperparameters are found by maximizing the Log Marginal Likelihood function. So, they are not guessed, but these are the hyperparameters that best describe the data. So finding the best value is basically the optimization process?
Choosing Non-Optimal Hyperparameters: Why would someone use the version without the best hyperparameters? One reason I’ve encountered is that optimization can take too long with the same data. Are there other reasons to choose non-optimizable models?
Best Practices: Are there standard practices or guidelines for when to use each version of the GPR models? Or one usually starts with non-optimizable version and then go to the optimizable one if the results are not sufficient enough?
Output Differences: Both versions output forecasted data and the covariance matrix. Does the optimized version provide any additional information or benefits?
Plotting Covariance: How can I plot the covariance that quantifies the uncertainty of my predictions? This is one of the main advantages of using GPR, and I want to visualize it. I have found some sample code here https://ch.mathworks.com/help/stats/gaussian-process-regression-models.html But I am a bit confused, because my input is a table actually, so I do not have only one predictor but a table. So far I have always used “predictedData = trainedModel.predictFcn(Table); but this gives me only the forecasted data without the 95% prediction intervals
I’m generally confused about the differences between both versions and the results they produce. Any insights or explanations would be greatly appreciated.
Thanks in advance!
Best regards,
Georgi regression, optimization MATLAB Answers — New Questions
compute sum of polynomial square
Can I solve sum of polynomial square in a one variable ?
e.g. (1+x)^2 + (3+x)^2 +2*x =
Can Mathlab compute result as x variable ?Can I solve sum of polynomial square in a one variable ?
e.g. (1+x)^2 + (3+x)^2 +2*x =
Can Mathlab compute result as x variable ? Can I solve sum of polynomial square in a one variable ?
e.g. (1+x)^2 + (3+x)^2 +2*x =
Can Mathlab compute result as x variable ? polynomial square MATLAB Answers — New Questions
Wrong font size for Default Axes Titles with set(groot,”DefaultAxesTitleFontSize”)
Hello,
I am trying to define some predefined figure styles by changing groot properties (I was inspired by the Default Property Values documentation). I am using the following syntax:
set(groot,"DefaultObjectTypePropertyName")
This works well for many cases, but I am unable to set correctly the Default Axes Title Font Size when I use LaTeX fonts.
I am using the "mwa_cmr10" font found in <MATLAB root folder>R2023asysfontsttf.
I actually can set the font size to the value I want, but after I divide it by 10.
Here is a MWE.
%% Dummy data
x1 = linspace(0,1,100);
y1 = sin(2*pi*x1);
x2 = linspace(0,1,100);
y2 = sin(4*pi*x2);
set( groot, "DefaultAxesTitleFontSize", 16/10 ); % Dividing the font size I want (16) by 10
set( groot, "DefaultTextInterpreter", "latex" );
set( groot, "DefaultAxesFontName", "mwa_cmr10" );
figure;
hold all
plot(x1, y1);
plot(x2, y2);
xlabel("X")
ylabel("Y")
title("Title")
legend;
gca().Title.FontSize % The font size is 16, but should be 1.6 as it is 16/10
I do not know if this is a bug or not, maybe it has something to do with the 10 in the "mwa_cmr10" font?
I wanted to ask this here first before signaling it as a bug. I think some more obscure graphics properties can be tricky to handle.
Thank you for your attention.
Regards,
Pedro
Edit: stylingHello,
I am trying to define some predefined figure styles by changing groot properties (I was inspired by the Default Property Values documentation). I am using the following syntax:
set(groot,"DefaultObjectTypePropertyName")
This works well for many cases, but I am unable to set correctly the Default Axes Title Font Size when I use LaTeX fonts.
I am using the "mwa_cmr10" font found in <MATLAB root folder>R2023asysfontsttf.
I actually can set the font size to the value I want, but after I divide it by 10.
Here is a MWE.
%% Dummy data
x1 = linspace(0,1,100);
y1 = sin(2*pi*x1);
x2 = linspace(0,1,100);
y2 = sin(4*pi*x2);
set( groot, "DefaultAxesTitleFontSize", 16/10 ); % Dividing the font size I want (16) by 10
set( groot, "DefaultTextInterpreter", "latex" );
set( groot, "DefaultAxesFontName", "mwa_cmr10" );
figure;
hold all
plot(x1, y1);
plot(x2, y2);
xlabel("X")
ylabel("Y")
title("Title")
legend;
gca().Title.FontSize % The font size is 16, but should be 1.6 as it is 16/10
I do not know if this is a bug or not, maybe it has something to do with the 10 in the "mwa_cmr10" font?
I wanted to ask this here first before signaling it as a bug. I think some more obscure graphics properties can be tricky to handle.
Thank you for your attention.
Regards,
Pedro
Edit: styling Hello,
I am trying to define some predefined figure styles by changing groot properties (I was inspired by the Default Property Values documentation). I am using the following syntax:
set(groot,"DefaultObjectTypePropertyName")
This works well for many cases, but I am unable to set correctly the Default Axes Title Font Size when I use LaTeX fonts.
I am using the "mwa_cmr10" font found in <MATLAB root folder>R2023asysfontsttf.
I actually can set the font size to the value I want, but after I divide it by 10.
Here is a MWE.
%% Dummy data
x1 = linspace(0,1,100);
y1 = sin(2*pi*x1);
x2 = linspace(0,1,100);
y2 = sin(4*pi*x2);
set( groot, "DefaultAxesTitleFontSize", 16/10 ); % Dividing the font size I want (16) by 10
set( groot, "DefaultTextInterpreter", "latex" );
set( groot, "DefaultAxesFontName", "mwa_cmr10" );
figure;
hold all
plot(x1, y1);
plot(x2, y2);
xlabel("X")
ylabel("Y")
title("Title")
legend;
gca().Title.FontSize % The font size is 16, but should be 1.6 as it is 16/10
I do not know if this is a bug or not, maybe it has something to do with the 10 in the "mwa_cmr10" font?
I wanted to ask this here first before signaling it as a bug. I think some more obscure graphics properties can be tricky to handle.
Thank you for your attention.
Regards,
Pedro
Edit: styling graphics, latex MATLAB Answers — New Questions
[Stateflow : Request for help] The question about extra time step for transition
Hello. I’m a newbie to stateflow in Simulink. Around 2-3 days ago, I figured out someting making me confused a little bit. I set the sample time in Model Configuration Paremeters to 0.01 sec, as shown in Fig 1. Fig 2. to Fig 3. shows the working process at each time step of the system wroted by myself, which consists of 3 parellel parent states. All of them are related to time counters (timer1, timer2 and timer3, which are local variables).
Fig 1. Sample time setting in Configuration Parameters (0.01 sec)
Fig 2. T = 0.000 sec
Fig 3. T = 0.010 sec
Fig 4. T = 0.020 sec
Fig 5. T = 0.030 sec
Fig 6. T = 0.040 sec
*Note that T is executing time, shown at the lowest of each picture.
According to the Fig 5., when T = 0.030 sec and all values of timer counters (timer1, timer2 and timer3) already reached to "3", there were only the buttom parent state, "Sys3_condition_action", that moved itself from "Off" state to "On" state, although all transitions between "On" and "Off" were satisfied with their timer counts (equal to 3). While the upper and middle parent states, Sys1_during and Sys2_cyclic_updating respectively, would move to "On" state in the next time step, at T = 0.040 sec (shown by Fig 6.).
At the first, I thought that all parents state should work in the same manner (all transitions should occur at T = 0.030 sec). But in fact, it’s not what I understand. Could someone explain to me why they worked differently?
Best Regards.Hello. I’m a newbie to stateflow in Simulink. Around 2-3 days ago, I figured out someting making me confused a little bit. I set the sample time in Model Configuration Paremeters to 0.01 sec, as shown in Fig 1. Fig 2. to Fig 3. shows the working process at each time step of the system wroted by myself, which consists of 3 parellel parent states. All of them are related to time counters (timer1, timer2 and timer3, which are local variables).
Fig 1. Sample time setting in Configuration Parameters (0.01 sec)
Fig 2. T = 0.000 sec
Fig 3. T = 0.010 sec
Fig 4. T = 0.020 sec
Fig 5. T = 0.030 sec
Fig 6. T = 0.040 sec
*Note that T is executing time, shown at the lowest of each picture.
According to the Fig 5., when T = 0.030 sec and all values of timer counters (timer1, timer2 and timer3) already reached to "3", there were only the buttom parent state, "Sys3_condition_action", that moved itself from "Off" state to "On" state, although all transitions between "On" and "Off" were satisfied with their timer counts (equal to 3). While the upper and middle parent states, Sys1_during and Sys2_cyclic_updating respectively, would move to "On" state in the next time step, at T = 0.040 sec (shown by Fig 6.).
At the first, I thought that all parents state should work in the same manner (all transitions should occur at T = 0.030 sec). But in fact, it’s not what I understand. Could someone explain to me why they worked differently?
Best Regards. Hello. I’m a newbie to stateflow in Simulink. Around 2-3 days ago, I figured out someting making me confused a little bit. I set the sample time in Model Configuration Paremeters to 0.01 sec, as shown in Fig 1. Fig 2. to Fig 3. shows the working process at each time step of the system wroted by myself, which consists of 3 parellel parent states. All of them are related to time counters (timer1, timer2 and timer3, which are local variables).
Fig 1. Sample time setting in Configuration Parameters (0.01 sec)
Fig 2. T = 0.000 sec
Fig 3. T = 0.010 sec
Fig 4. T = 0.020 sec
Fig 5. T = 0.030 sec
Fig 6. T = 0.040 sec
*Note that T is executing time, shown at the lowest of each picture.
According to the Fig 5., when T = 0.030 sec and all values of timer counters (timer1, timer2 and timer3) already reached to "3", there were only the buttom parent state, "Sys3_condition_action", that moved itself from "Off" state to "On" state, although all transitions between "On" and "Off" were satisfied with their timer counts (equal to 3). While the upper and middle parent states, Sys1_during and Sys2_cyclic_updating respectively, would move to "On" state in the next time step, at T = 0.040 sec (shown by Fig 6.).
At the first, I thought that all parents state should work in the same manner (all transitions should occur at T = 0.030 sec). But in fact, it’s not what I understand. Could someone explain to me why they worked differently?
Best Regards. stateflow, simulink MATLAB Answers — New Questions
How to choose a single element randomly from a vector
A=[2 3 4 5];
How do I choose any one variable from the vector A randomly each time.A=[2 3 4 5];
How do I choose any one variable from the vector A randomly each time. A=[2 3 4 5];
How do I choose any one variable from the vector A randomly each time. rand, randi MATLAB Answers — New Questions
I want to export geometry generated in MATLAB, in abaqus for further FEM analysis. Can any one suggest me how to export geometry for abaqus ?
I had prepared a tool to generate user defined geometry of brick. Now I want to to use that tool to generate various geometries for analysis to be done in abaqus.I had prepared a tool to generate user defined geometry of brick. Now I want to to use that tool to generate various geometries for analysis to be done in abaqus. I had prepared a tool to generate user defined geometry of brick. Now I want to to use that tool to generate various geometries for analysis to be done in abaqus. geometry export, export to abaqus MATLAB Answers — New Questions
Unable to recreate k-means clustering example
As per title, i’m unable to recreate the following example (k-means clustering – MATLAB kmeans – MathWorks Italia), the k means clustering of Fisher’s iris dataset.
This is the code in question
load fisheriris
X = meas(:,3:4);
figure;
plot(X(:,1),X(:,2),’k*’,’MarkerSize’,5);
title ‘Fisher”s Iris Data’;
xlabel ‘Petal Lengths (cm)’;
ylabel ‘Petal Widths (cm)’
rng(1); % For reproducibility
[idx,C] = kmeans(X,3);
x1 = min(X(:,1)):0.01:max(X(:,1));
x2 = min(X(:,2)):0.01:max(X(:,2));
[x1G,x2G] = meshgrid(x1,x2);
XGrid = [x1G(:),x2G(:)]; % Defines a fine grid on the plot
idx2Region = kmeans(XGrid,3,’MaxIter’,1,’Start’,C);
figure;
gscatter(XGrid(:,1),XGrid(:,2),idx2Region,…
[0,0.75,0.75;0.75,0,0.75;0.75,0.75,0],’..’);
hold on;
plot(X(:,1),X(:,2),’k*’,’MarkerSize’,5);
title ‘Fisher”s Iris Data’;
xlabel ‘Petal Lengths (cm)’;
ylabel ‘Petal Widths (cm)’;
legend(‘Region 1′,’Region 2′,’Region 3′,’Data’,’Location’,’SouthEast’);
hold off;
An error occurs at line 16 involving the kmeans function
Error in Kmeans_dimostrativo (line 16)
[idx,C] = kmeans(X,3);
So i started digging in the kmeans function
function [label, mu, energy] = kmeans(X, m)
% Perform kmeans clustering.
% Input:
% X: d x n data matrix
% m: initialization parameter
% Output:
% label: 1 x n sample labels
% mu: d x k center of clusters
% energy: optimization target value
% Written by Mo Chen (sth4nth@gmail.com).
label = init(X, m);
n = numel(label);
idx = 1:n;
last = zeros(1,n);
while any(label ~= last)
[~,~,last(:)] = unique(label); % remove empty clusters
mu = X*normalize(sparse(idx,last,1),1); % compute cluster centers
[val,label] = min(dot(mu,mu,1)’/2-mu’*X,[],1); % assign sample labels
end
energy = dot(X(:),X(:),1)+2*sum(val);
function label = init(X, m)
[d,n] = size(X);
if numel(m) == 1 % random initialization
mu = X(:,randperm(n,m));
[~,label] = min(dot(mu,mu,1)’/2-mu’*X,[],1);
elseif all(size(m) == [1,n]) % init with labels
label = m;
elseif size(m,1) == d % init with seeds (centers)
[~,label] = min(dot(m,m,1)’/2-m’*X,[],1);
end
The following error messages are the ones relative to the kmeans function
Error using randperm
K must be less than or equal to N.
Error in kmeans>init (line 25)
mu = X(:,randperm(n,m));
Error in kmeans (line 11)
label = init(X, m);
I honestly don’t know the reason for the errors, i took the script directly from the website.As per title, i’m unable to recreate the following example (k-means clustering – MATLAB kmeans – MathWorks Italia), the k means clustering of Fisher’s iris dataset.
This is the code in question
load fisheriris
X = meas(:,3:4);
figure;
plot(X(:,1),X(:,2),’k*’,’MarkerSize’,5);
title ‘Fisher”s Iris Data’;
xlabel ‘Petal Lengths (cm)’;
ylabel ‘Petal Widths (cm)’
rng(1); % For reproducibility
[idx,C] = kmeans(X,3);
x1 = min(X(:,1)):0.01:max(X(:,1));
x2 = min(X(:,2)):0.01:max(X(:,2));
[x1G,x2G] = meshgrid(x1,x2);
XGrid = [x1G(:),x2G(:)]; % Defines a fine grid on the plot
idx2Region = kmeans(XGrid,3,’MaxIter’,1,’Start’,C);
figure;
gscatter(XGrid(:,1),XGrid(:,2),idx2Region,…
[0,0.75,0.75;0.75,0,0.75;0.75,0.75,0],’..’);
hold on;
plot(X(:,1),X(:,2),’k*’,’MarkerSize’,5);
title ‘Fisher”s Iris Data’;
xlabel ‘Petal Lengths (cm)’;
ylabel ‘Petal Widths (cm)’;
legend(‘Region 1′,’Region 2′,’Region 3′,’Data’,’Location’,’SouthEast’);
hold off;
An error occurs at line 16 involving the kmeans function
Error in Kmeans_dimostrativo (line 16)
[idx,C] = kmeans(X,3);
So i started digging in the kmeans function
function [label, mu, energy] = kmeans(X, m)
% Perform kmeans clustering.
% Input:
% X: d x n data matrix
% m: initialization parameter
% Output:
% label: 1 x n sample labels
% mu: d x k center of clusters
% energy: optimization target value
% Written by Mo Chen (sth4nth@gmail.com).
label = init(X, m);
n = numel(label);
idx = 1:n;
last = zeros(1,n);
while any(label ~= last)
[~,~,last(:)] = unique(label); % remove empty clusters
mu = X*normalize(sparse(idx,last,1),1); % compute cluster centers
[val,label] = min(dot(mu,mu,1)’/2-mu’*X,[],1); % assign sample labels
end
energy = dot(X(:),X(:),1)+2*sum(val);
function label = init(X, m)
[d,n] = size(X);
if numel(m) == 1 % random initialization
mu = X(:,randperm(n,m));
[~,label] = min(dot(mu,mu,1)’/2-mu’*X,[],1);
elseif all(size(m) == [1,n]) % init with labels
label = m;
elseif size(m,1) == d % init with seeds (centers)
[~,label] = min(dot(m,m,1)’/2-m’*X,[],1);
end
The following error messages are the ones relative to the kmeans function
Error using randperm
K must be less than or equal to N.
Error in kmeans>init (line 25)
mu = X(:,randperm(n,m));
Error in kmeans (line 11)
label = init(X, m);
I honestly don’t know the reason for the errors, i took the script directly from the website. As per title, i’m unable to recreate the following example (k-means clustering – MATLAB kmeans – MathWorks Italia), the k means clustering of Fisher’s iris dataset.
This is the code in question
load fisheriris
X = meas(:,3:4);
figure;
plot(X(:,1),X(:,2),’k*’,’MarkerSize’,5);
title ‘Fisher”s Iris Data’;
xlabel ‘Petal Lengths (cm)’;
ylabel ‘Petal Widths (cm)’
rng(1); % For reproducibility
[idx,C] = kmeans(X,3);
x1 = min(X(:,1)):0.01:max(X(:,1));
x2 = min(X(:,2)):0.01:max(X(:,2));
[x1G,x2G] = meshgrid(x1,x2);
XGrid = [x1G(:),x2G(:)]; % Defines a fine grid on the plot
idx2Region = kmeans(XGrid,3,’MaxIter’,1,’Start’,C);
figure;
gscatter(XGrid(:,1),XGrid(:,2),idx2Region,…
[0,0.75,0.75;0.75,0,0.75;0.75,0.75,0],’..’);
hold on;
plot(X(:,1),X(:,2),’k*’,’MarkerSize’,5);
title ‘Fisher”s Iris Data’;
xlabel ‘Petal Lengths (cm)’;
ylabel ‘Petal Widths (cm)’;
legend(‘Region 1′,’Region 2′,’Region 3′,’Data’,’Location’,’SouthEast’);
hold off;
An error occurs at line 16 involving the kmeans function
Error in Kmeans_dimostrativo (line 16)
[idx,C] = kmeans(X,3);
So i started digging in the kmeans function
function [label, mu, energy] = kmeans(X, m)
% Perform kmeans clustering.
% Input:
% X: d x n data matrix
% m: initialization parameter
% Output:
% label: 1 x n sample labels
% mu: d x k center of clusters
% energy: optimization target value
% Written by Mo Chen (sth4nth@gmail.com).
label = init(X, m);
n = numel(label);
idx = 1:n;
last = zeros(1,n);
while any(label ~= last)
[~,~,last(:)] = unique(label); % remove empty clusters
mu = X*normalize(sparse(idx,last,1),1); % compute cluster centers
[val,label] = min(dot(mu,mu,1)’/2-mu’*X,[],1); % assign sample labels
end
energy = dot(X(:),X(:),1)+2*sum(val);
function label = init(X, m)
[d,n] = size(X);
if numel(m) == 1 % random initialization
mu = X(:,randperm(n,m));
[~,label] = min(dot(mu,mu,1)’/2-mu’*X,[],1);
elseif all(size(m) == [1,n]) % init with labels
label = m;
elseif size(m,1) == d % init with seeds (centers)
[~,label] = min(dot(m,m,1)’/2-m’*X,[],1);
end
The following error messages are the ones relative to the kmeans function
Error using randperm
K must be less than or equal to N.
Error in kmeans>init (line 25)
mu = X(:,randperm(n,m));
Error in kmeans (line 11)
label = init(X, m);
I honestly don’t know the reason for the errors, i took the script directly from the website. k-means clustering, statistics MATLAB Answers — New Questions
How to receive an full binary data using mqtt callback function
I have sent a binary blob (for example 4k) from C++ app by mqtt mosquitto. In matlab only a chunk of sending data is received. This chunk is limited by the first zero byte which is in the blob. How to receive a full message?
Short code snippest below:
mqttClient = mqttclient("tcp://127.0.0.1");
mySub = subscribe(mqttClient, "topic", Callback=@MsgCallvBack)
%
function MsgCallvBack(topic, data)
fprintf("topic=%s, data size=%un", topic, numel(char(data)));
endI have sent a binary blob (for example 4k) from C++ app by mqtt mosquitto. In matlab only a chunk of sending data is received. This chunk is limited by the first zero byte which is in the blob. How to receive a full message?
Short code snippest below:
mqttClient = mqttclient("tcp://127.0.0.1");
mySub = subscribe(mqttClient, "topic", Callback=@MsgCallvBack)
%
function MsgCallvBack(topic, data)
fprintf("topic=%s, data size=%un", topic, numel(char(data)));
end I have sent a binary blob (for example 4k) from C++ app by mqtt mosquitto. In matlab only a chunk of sending data is received. This chunk is limited by the first zero byte which is in the blob. How to receive a full message?
Short code snippest below:
mqttClient = mqttclient("tcp://127.0.0.1");
mySub = subscribe(mqttClient, "topic", Callback=@MsgCallvBack)
%
function MsgCallvBack(topic, data)
fprintf("topic=%s, data size=%un", topic, numel(char(data)));
end mqtt MATLAB Answers — New Questions
Error using fit>iFit Too many start points. You need 2 start points for this model
Dear all
I have created this function to do implement a fitting to some data:
function [fitresult,gof]=Fit(s,mag,in_mag,in_cond)
if round(in_mag(1))==1 && round(in_mag(end))==-1
fitting_function=’cos(2*atan(exp((x-x0)/Delta)))’;
elseif round(in_mag(1))==-1 && round(in_mag(end))==1
fitting_function=’cos(2*atan(exp(-(x-x0)/Delta)))’;
end
ft=fittype(fitting_function,’independent’,’x’,’dependent’,’y’);
opts=fitoptions(‘Method’,’NonlinearLeastSquares’);
opts.Display=’Off’;
opts.StartPoint=in_cond;
[fitresult,gof]=fit(space,mag,ft,opts);
When trying to fit some data:
initial_conditions=[1 1];
[fitresult,gof]=Fit(some_column_vector,other_column_vector,another_column_vector,initial_conditions’);
where "some_column_vector", "other_column_vector", and "another_column_vector" are 1000×1 arrays.
I received the following error:
Error using fit>iFit
Too many start points. You need 2 start points for this model.
Any ideas on what could be happening here?Dear all
I have created this function to do implement a fitting to some data:
function [fitresult,gof]=Fit(s,mag,in_mag,in_cond)
if round(in_mag(1))==1 && round(in_mag(end))==-1
fitting_function=’cos(2*atan(exp((x-x0)/Delta)))’;
elseif round(in_mag(1))==-1 && round(in_mag(end))==1
fitting_function=’cos(2*atan(exp(-(x-x0)/Delta)))’;
end
ft=fittype(fitting_function,’independent’,’x’,’dependent’,’y’);
opts=fitoptions(‘Method’,’NonlinearLeastSquares’);
opts.Display=’Off’;
opts.StartPoint=in_cond;
[fitresult,gof]=fit(space,mag,ft,opts);
When trying to fit some data:
initial_conditions=[1 1];
[fitresult,gof]=Fit(some_column_vector,other_column_vector,another_column_vector,initial_conditions’);
where "some_column_vector", "other_column_vector", and "another_column_vector" are 1000×1 arrays.
I received the following error:
Error using fit>iFit
Too many start points. You need 2 start points for this model.
Any ideas on what could be happening here? Dear all
I have created this function to do implement a fitting to some data:
function [fitresult,gof]=Fit(s,mag,in_mag,in_cond)
if round(in_mag(1))==1 && round(in_mag(end))==-1
fitting_function=’cos(2*atan(exp((x-x0)/Delta)))’;
elseif round(in_mag(1))==-1 && round(in_mag(end))==1
fitting_function=’cos(2*atan(exp(-(x-x0)/Delta)))’;
end
ft=fittype(fitting_function,’independent’,’x’,’dependent’,’y’);
opts=fitoptions(‘Method’,’NonlinearLeastSquares’);
opts.Display=’Off’;
opts.StartPoint=in_cond;
[fitresult,gof]=fit(space,mag,ft,opts);
When trying to fit some data:
initial_conditions=[1 1];
[fitresult,gof]=Fit(some_column_vector,other_column_vector,another_column_vector,initial_conditions’);
where "some_column_vector", "other_column_vector", and "another_column_vector" are 1000×1 arrays.
I received the following error:
Error using fit>iFit
Too many start points. You need 2 start points for this model.
Any ideas on what could be happening here? fitting MATLAB Answers — New Questions
How to approximate a measured frequency response function to the sum of 2 1DOF vibration FRF functions?
I have measured the vibration frequency response function a structure (blue curves in the following diagrams).
By using command
SYSD=tfest(SYSFR,5);
% SYSFR is the measured result (imported beginning from 25 Hz, as the
% low-frequency area is not accurate). SYSD is the approximated FRF.
I can approximate the measured FRF to a transfer function with 5 poles (the least amount of poles needed for the approximation to be accurate enough).
Since the measured FRF curves clearly have the shape of 2DOF vibration system, I want to further approximate the transfer function as the sum of two 1DOF vibration transfer functions, and draw both 1DOF vibration FRF curves in diagram (like in the following diagram). How to achieve this? The approximation can either be done based on measured FRF or approximated 5-pole FRF.I have measured the vibration frequency response function a structure (blue curves in the following diagrams).
By using command
SYSD=tfest(SYSFR,5);
% SYSFR is the measured result (imported beginning from 25 Hz, as the
% low-frequency area is not accurate). SYSD is the approximated FRF.
I can approximate the measured FRF to a transfer function with 5 poles (the least amount of poles needed for the approximation to be accurate enough).
Since the measured FRF curves clearly have the shape of 2DOF vibration system, I want to further approximate the transfer function as the sum of two 1DOF vibration transfer functions, and draw both 1DOF vibration FRF curves in diagram (like in the following diagram). How to achieve this? The approximation can either be done based on measured FRF or approximated 5-pole FRF. I have measured the vibration frequency response function a structure (blue curves in the following diagrams).
By using command
SYSD=tfest(SYSFR,5);
% SYSFR is the measured result (imported beginning from 25 Hz, as the
% low-frequency area is not accurate). SYSD is the approximated FRF.
I can approximate the measured FRF to a transfer function with 5 poles (the least amount of poles needed for the approximation to be accurate enough).
Since the measured FRF curves clearly have the shape of 2DOF vibration system, I want to further approximate the transfer function as the sum of two 1DOF vibration transfer functions, and draw both 1DOF vibration FRF curves in diagram (like in the following diagram). How to achieve this? The approximation can either be done based on measured FRF or approximated 5-pole FRF. frequency, transfer function MATLAB Answers — New Questions
I want my code conditions checked every 1/2000th of a second, how can I do that?
I have multiple if conditions that is being checked to provide pulse and gate signals to switches I want that these conditions be checked every 1/2000th of a second. I tried for loop and tic toc methods but the condition is being checked at a frequency of 200kHz. I am not able to solve this issue kindly helpI have multiple if conditions that is being checked to provide pulse and gate signals to switches I want that these conditions be checked every 1/2000th of a second. I tried for loop and tic toc methods but the condition is being checked at a frequency of 200kHz. I am not able to solve this issue kindly help I have multiple if conditions that is being checked to provide pulse and gate signals to switches I want that these conditions be checked every 1/2000th of a second. I tried for loop and tic toc methods but the condition is being checked at a frequency of 200kHz. I am not able to solve this issue kindly help loop condition, if statement, for loop MATLAB Answers — New Questions