Hello and sorry for this long query. I am not very good at this , I have to admit it.
I am encountering an issue when trying to apply custom boundary conditions in a PDE model using MATLAB. Specifically, I have a custom boundary condition function mybc1 that works correctly in another code when using structuralBC to apply displacement conditions to vertices. However, in my current code, when I try to use applyBoundaryCondition to set mixed boundary conditions on specific faces, it seems that the solver does not enter the custom boundary condition function mybc1.
Here’s a brief description of what I am trying to do:
Model Definition: I have defined a PDE model with certain geometry and mesh.
Boundary Condition Application: I attempt to apply mixed boundary conditions on specific faces using applyBoundaryCondition with the u parameter set to a custom function handle @(location, state)mybc1(location, state, dis, ts, in).
Issue: It appears that the solver does not call the mybc1 function during the solution process. I verified this by setting breakpoints inside mybc1, which are never hit.
The mybc1 function interpolates the acceleration data (dis) at the times specified by the solver (state.time) and returns the corresponding values. This function has been tested and works correctly in another program with structuralBC.

code as follows:

%____________Define dimensions of the plate________________________________
len = 1.22; % length in x direction
width = 1.22; % Width in y direction
sq_side = 0.10; % Side length of the squares in corners
% Create the geometry description matrix___________________________________
% —–Define the outer boundary——————————————-
outer_boundary = [3, 4, 0, len, len, 0, 0, 0, width, width]’;
%——Define squares in the corners—————————————
square_1 = [3, 4, 0, sq_side, sq_side, 0, 0, 0, sq_side, sq_side]’;
square_2 = [3, 4, len – sq_side, len, len, len – sq_side, 0, 0, sq_side, sq_side]’;
square_3 = [3, 4, len – sq_side, len, len, len – sq_side, width – sq_side, width – sq_side, width, width]’;
square_4 = [3, 4, 0, sq_side, sq_side, 0, width – sq_side, width – sq_side, width, width]’;
%——Combine all the geometries——————————————
gdm = [outer_boundary, square_1, square_2, square_3, square_4];
%——Define the names for each region————————————
ns = (char(‘R1′,’R2′,’R3′,’R4′,’R5′))’;
sf = ‘R1+R2+R3+R4+R5’;
%——-Create the geometry————————————————
model = createpde(2);
geometryFromEdges(model, decsg(gdm, sf, ns));
%——-Generate the mesh and plot the geometry—————————-
msh = generateMesh(model, ‘Hmax’, 0.1); % Use ‘Hmax’ to control mesh density
figure
pdemesh(model)
%_____________Create the different regions of the model___________________
% Obtain nodes and elements of the mesh
[p,e,t] = meshToPet(model.Mesh);
% ———-Plot with Faces labels—————————————-
figure;
pdegplot(model, ‘FaceLabels’, ‘on’, ‘FaceAlpha’, 0.5);
figure
pdemesh(model, ‘NodeLabels’, ‘on’);
% Definition of the constants
E = 4E9;
h_thick = 0.05;
nu = 0.3;
mass = 100;
D = E*h_thick^3/(12*(1-nu)^2);
% Now we create the PDE systems as symbolic equations_____________________
syms pres
syms u1(x,y,t) u2(x,y,t)
pdeeq = [-laplacian(u1,[x y])+u2; D*laplacian(u2,[x y])+ mass*diff(u1,t,t)-pres];
symcoeffs = pdeCoefficients(pdeeq,[u1,u2],’Symbolic’,true);
c2=symcoeffs.c;
m2=symcoeffs.m;
%——Display the symbolic coefficients———————————–
structfun(@disp,symcoeffs);
symcoeffs=subs(symcoeffs,pres,1);
coeffs=pdeCoefficientsToDouble(symcoeffs);
%——pass these coefficients to the pde model—————————-
specifyCoefficients(model,’m’,coeffs.m,’d’,coeffs.d,’c’,…
coeffs.c,’a’,coeffs.a,’f’,coeffs.f);
% INITIAL CONDITIONS ——————————————————
setInitialConditions(model,[0;0],[0;0]);
%USE THE DATA FROM ASHMOLEAN AS BOUNDARY CONDITIONS————————
%load the vector containing the acceleration values
load(totaldata.mat’, ‘DATA’,…
‘TIME’,’fs’);
ac=DATA;
% ——-Create vector time, not starting in 0
num_samples = size(DATA, 1);
dt = 1 / fs; % interval between measurements
t = (1:num_samples) * dt; % vector of times
% TRANSFORM ACCELERATIONS INTO DISPLACEMENTS_______________________________
acvedi=AccVelDis(962.54,1005.11,992.35,ac,fs,t);
dis=acvedi{3};
t=[0,t]; % I start the time at 0
% Boundary conditions in the faces using the displacement vectors
in=0; %counter variable
Faces=[1,2,3,5]; %Vector with faces numbers
for jk = 1:numel(Faces)
in=in+1;%add 1 in every loop to access the different measurmentes contained in the matrix acceleration
face= Faces(jk);
dis1=dis(:,in);

applyBoundaryCondition(model,"mixed","Face",[face,face],"u",@(location,state)mybc1(location,state,dis1,t),"EquationIndex",1,"q",[0 0],"g",0);
%applyBoundaryCondition(model,"dirichlet","Face",face,"h",[0 0 ; 0 1],"r",bcfunc)
end
in=0;
tim=[t(2) t(3)]; %tim represents the time of interest to solve the pde
res=solvepde(model, tim);
% Access the solution at the nodal locations
sol=res.NodalSolution;

And this is the function that works in another code:

function bcMatrix = mybc1(location,state,acs,ti)
%This function is use to extrapolate the values of acceleration for the
%time instances chosen by the system to do the integration of the system.
%as imput it requires the measurements matrix and the corresponding time
%measurements. As output, it provices the extrapoleted values for the state
%time

T=state.time;
% Check if T is NaN and assign the previous value if true
if isnan(T)
vq = NaN(size(location.x)); %this is what documentation says????
else
ac = acs; %gets the value of the corresponding acceleration
vq = interp1(ti, ac, T); %interpolates it
end
bcMatrix = vq;
end

I would really appreciate any help/guidanceHello and sorry for this long query. I am not very good at this , I have to admit it.
I am encountering an issue when trying to apply custom boundary conditions in a PDE model using MATLAB. Specifically, I have a custom boundary condition function mybc1 that works correctly in another code when using structuralBC to apply displacement conditions to vertices. However, in my current code, when I try to use applyBoundaryCondition to set mixed boundary conditions on specific faces, it seems that the solver does not enter the custom boundary condition function mybc1.
Here’s a brief description of what I am trying to do:
Model Definition: I have defined a PDE model with certain geometry and mesh.
Boundary Condition Application: I attempt to apply mixed boundary conditions on specific faces using applyBoundaryCondition with the u parameter set to a custom function handle @(location, state)mybc1(location, state, dis, ts, in).
Issue: It appears that the solver does not call the mybc1 function during the solution process. I verified this by setting breakpoints inside mybc1, which are never hit.
The mybc1 function interpolates the acceleration data (dis) at the times specified by the solver (state.time) and returns the corresponding values. This function has been tested and works correctly in another program with structuralBC.

code as follows:

%____________Define dimensions of the plate________________________________
len = 1.22; % length in x direction
width = 1.22; % Width in y direction
sq_side = 0.10; % Side length of the squares in corners
% Create the geometry description matrix___________________________________
% —–Define the outer boundary——————————————-
outer_boundary = [3, 4, 0, len, len, 0, 0, 0, width, width]’;
%——Define squares in the corners—————————————
square_1 = [3, 4, 0, sq_side, sq_side, 0, 0, 0, sq_side, sq_side]’;
square_2 = [3, 4, len – sq_side, len, len, len – sq_side, 0, 0, sq_side, sq_side]’;
square_3 = [3, 4, len – sq_side, len, len, len – sq_side, width – sq_side, width – sq_side, width, width]’;
square_4 = [3, 4, 0, sq_side, sq_side, 0, width – sq_side, width – sq_side, width, width]’;
%——Combine all the geometries——————————————
gdm = [outer_boundary, square_1, square_2, square_3, square_4];
%——Define the names for each region————————————
ns = (char(‘R1′,’R2′,’R3′,’R4′,’R5′))’;
sf = ‘R1+R2+R3+R4+R5’;
%——-Create the geometry————————————————
model = createpde(2);
geometryFromEdges(model, decsg(gdm, sf, ns));
%——-Generate the mesh and plot the geometry—————————-
msh = generateMesh(model, ‘Hmax’, 0.1); % Use ‘Hmax’ to control mesh density
figure
pdemesh(model)
%_____________Create the different regions of the model___________________
% Obtain nodes and elements of the mesh
[p,e,t] = meshToPet(model.Mesh);
% ———-Plot with Faces labels—————————————-
figure;
pdegplot(model, ‘FaceLabels’, ‘on’, ‘FaceAlpha’, 0.5);
figure
pdemesh(model, ‘NodeLabels’, ‘on’);
% Definition of the constants
E = 4E9;
h_thick = 0.05;
nu = 0.3;
mass = 100;
D = E*h_thick^3/(12*(1-nu)^2);
% Now we create the PDE systems as symbolic equations_____________________
syms pres
syms u1(x,y,t) u2(x,y,t)
pdeeq = [-laplacian(u1,[x y])+u2; D*laplacian(u2,[x y])+ mass*diff(u1,t,t)-pres];
symcoeffs = pdeCoefficients(pdeeq,[u1,u2],’Symbolic’,true);
c2=symcoeffs.c;
m2=symcoeffs.m;
%——Display the symbolic coefficients———————————–
structfun(@disp,symcoeffs);
symcoeffs=subs(symcoeffs,pres,1);
coeffs=pdeCoefficientsToDouble(symcoeffs);
%——pass these coefficients to the pde model—————————-
specifyCoefficients(model,’m’,coeffs.m,’d’,coeffs.d,’c’,…
coeffs.c,’a’,coeffs.a,’f’,coeffs.f);
% INITIAL CONDITIONS ——————————————————
setInitialConditions(model,[0;0],[0;0]);
%USE THE DATA FROM ASHMOLEAN AS BOUNDARY CONDITIONS————————
%load the vector containing the acceleration values
load(totaldata.mat’, ‘DATA’,…
‘TIME’,’fs’);
ac=DATA;
% ——-Create vector time, not starting in 0
num_samples = size(DATA, 1);
dt = 1 / fs; % interval between measurements
t = (1:num_samples) * dt; % vector of times
% TRANSFORM ACCELERATIONS INTO DISPLACEMENTS_______________________________
acvedi=AccVelDis(962.54,1005.11,992.35,ac,fs,t);
dis=acvedi{3};
t=[0,t]; % I start the time at 0
% Boundary conditions in the faces using the displacement vectors
in=0; %counter variable
Faces=[1,2,3,5]; %Vector with faces numbers
for jk = 1:numel(Faces)
in=in+1;%add 1 in every loop to access the different measurmentes contained in the matrix acceleration
face= Faces(jk);
dis1=dis(:,in);

applyBoundaryCondition(model,"mixed","Face",[face,face],"u",@(location,state)mybc1(location,state,dis1,t),"EquationIndex",1,"q",[0 0],"g",0);
%applyBoundaryCondition(model,"dirichlet","Face",face,"h",[0 0 ; 0 1],"r",bcfunc)
end
in=0;
tim=[t(2) t(3)]; %tim represents the time of interest to solve the pde
res=solvepde(model, tim);
% Access the solution at the nodal locations
sol=res.NodalSolution;

And this is the function that works in another code:

function bcMatrix = mybc1(location,state,acs,ti)
%This function is use to extrapolate the values of acceleration for the
%time instances chosen by the system to do the integration of the system.
%as imput it requires the measurements matrix and the corresponding time
%measurements. As output, it provices the extrapoleted values for the state
%time

T=state.time;
% Check if T is NaN and assign the previous value if true
if isnan(T)
vq = NaN(size(location.x)); %this is what documentation says????
else
ac = acs; %gets the value of the corresponding acceleration
vq = interp1(ti, ac, T); %interpolates it
end
bcMatrix = vq;
end

I would really appreciate any help/guidance Hello and sorry for this long query. I am not very good at this , I have to admit it.
I am encountering an issue when trying to apply custom boundary conditions in a PDE model using MATLAB. Specifically, I have a custom boundary condition function mybc1 that works correctly in another code when using structuralBC to apply displacement conditions to vertices. However, in my current code, when I try to use applyBoundaryCondition to set mixed boundary conditions on specific faces, it seems that the solver does not enter the custom boundary condition function mybc1.
Here’s a brief description of what I am trying to do:
Model Definition: I have defined a PDE model with certain geometry and mesh.
Boundary Condition Application: I attempt to apply mixed boundary conditions on specific faces using applyBoundaryCondition with the u parameter set to a custom function handle @(location, state)mybc1(location, state, dis, ts, in).
Issue: It appears that the solver does not call the mybc1 function during the solution process. I verified this by setting breakpoints inside mybc1, which are never hit.
The mybc1 function interpolates the acceleration data (dis) at the times specified by the solver (state.time) and returns the corresponding values. This function has been tested and works correctly in another program with structuralBC.

code as follows:

%____________Define dimensions of the plate________________________________
len = 1.22; % length in x direction
width = 1.22; % Width in y direction
sq_side = 0.10; % Side length of the squares in corners
% Create the geometry description matrix___________________________________
% —–Define the outer boundary——————————————-
outer_boundary = [3, 4, 0, len, len, 0, 0, 0, width, width]’;
%——Define squares in the corners—————————————
square_1 = [3, 4, 0, sq_side, sq_side, 0, 0, 0, sq_side, sq_side]’;
square_2 = [3, 4, len – sq_side, len, len, len – sq_side, 0, 0, sq_side, sq_side]’;
square_3 = [3, 4, len – sq_side, len, len, len – sq_side, width – sq_side, width – sq_side, width, width]’;
square_4 = [3, 4, 0, sq_side, sq_side, 0, width – sq_side, width – sq_side, width, width]’;
%——Combine all the geometries——————————————
gdm = [outer_boundary, square_1, square_2, square_3, square_4];
%——Define the names for each region————————————
ns = (char(‘R1′,’R2′,’R3′,’R4′,’R5′))’;
sf = ‘R1+R2+R3+R4+R5’;
%——-Create the geometry————————————————
model = createpde(2);
geometryFromEdges(model, decsg(gdm, sf, ns));
%——-Generate the mesh and plot the geometry—————————-
msh = generateMesh(model, ‘Hmax’, 0.1); % Use ‘Hmax’ to control mesh density
figure
pdemesh(model)
%_____________Create the different regions of the model___________________
% Obtain nodes and elements of the mesh
[p,e,t] = meshToPet(model.Mesh);
% ———-Plot with Faces labels—————————————-
figure;
pdegplot(model, ‘FaceLabels’, ‘on’, ‘FaceAlpha’, 0.5);
figure
pdemesh(model, ‘NodeLabels’, ‘on’);
% Definition of the constants
E = 4E9;
h_thick = 0.05;
nu = 0.3;
mass = 100;
D = E*h_thick^3/(12*(1-nu)^2);
% Now we create the PDE systems as symbolic equations_____________________
syms pres
syms u1(x,y,t) u2(x,y,t)
pdeeq = [-laplacian(u1,[x y])+u2; D*laplacian(u2,[x y])+ mass*diff(u1,t,t)-pres];
symcoeffs = pdeCoefficients(pdeeq,[u1,u2],’Symbolic’,true);
c2=symcoeffs.c;
m2=symcoeffs.m;
%——Display the symbolic coefficients———————————–
structfun(@disp,symcoeffs);
symcoeffs=subs(symcoeffs,pres,1);
coeffs=pdeCoefficientsToDouble(symcoeffs);
%——pass these coefficients to the pde model—————————-
specifyCoefficients(model,’m’,coeffs.m,’d’,coeffs.d,’c’,…
coeffs.c,’a’,coeffs.a,’f’,coeffs.f);
% INITIAL CONDITIONS ——————————————————
setInitialConditions(model,[0;0],[0;0]);
%USE THE DATA FROM ASHMOLEAN AS BOUNDARY CONDITIONS————————
%load the vector containing the acceleration values
load(totaldata.mat’, ‘DATA’,…
‘TIME’,’fs’);
ac=DATA;
% ——-Create vector time, not starting in 0
num_samples = size(DATA, 1);
dt = 1 / fs; % interval between measurements
t = (1:num_samples) * dt; % vector of times
% TRANSFORM ACCELERATIONS INTO DISPLACEMENTS_______________________________
acvedi=AccVelDis(962.54,1005.11,992.35,ac,fs,t);
dis=acvedi{3};
t=[0,t]; % I start the time at 0
% Boundary conditions in the faces using the displacement vectors
in=0; %counter variable
Faces=[1,2,3,5]; %Vector with faces numbers
for jk = 1:numel(Faces)
in=in+1;%add 1 in every loop to access the different measurmentes contained in the matrix acceleration
face= Faces(jk);
dis1=dis(:,in);

applyBoundaryCondition(model,"mixed","Face",[face,face],"u",@(location,state)mybc1(location,state,dis1,t),"EquationIndex",1,"q",[0 0],"g",0);
%applyBoundaryCondition(model,"dirichlet","Face",face,"h",[0 0 ; 0 1],"r",bcfunc)
end
in=0;
tim=[t(2) t(3)]; %tim represents the time of interest to solve the pde
res=solvepde(model, tim);
% Access the solution at the nodal locations
sol=res.NodalSolution;

And this is the function that works in another code:

function bcMatrix = mybc1(location,state,acs,ti)
%This function is use to extrapolate the values of acceleration for the
%time instances chosen by the system to do the integration of the system.
%as imput it requires the measurements matrix and the corresponding time
%measurements. As output, it provices the extrapoleted values for the state
%time

T=state.time;
% Check if T is NaN and assign the previous value if true
if isnan(T)
vq = NaN(size(location.x)); %this is what documentation says????
else
ac = acs; %gets the value of the corresponding acceleration
vq = interp1(ti, ac, T); %interpolates it
end
bcMatrix = vq;
end

I would really appreciate any help/guidance handles, pde, applyboundaryconditions, function, differential equations MATLAB Answers — New Questions

​

Hello, I have two matrices, M1 ans M2. They both have rows and columns fom 0 to 9 and a to f (16 * 16) matrix.
Mi matrix is filled with random numbers form 0 to 255 in hex decimal notation.

The second matrix M2 is empty but has the same row and column indexs. I need to fill M2 with the help of M1. For eample, in M2, at row 3 and colm 4 ( 34,the first entry of M1) will be 00 ( index of the first entry of M1). In M2 at row 6 col A, the entry will be 01, at row 2 col 9, the entry will be 02, and so on.
%% % Initialize M2 with zeros
M2 = zeros(16,16);
for row = 0:15
for col = 0:15
% get the value from M1
value= M1(row+1,col+1);

% Convert the values to row and col indices for M2

hexStr=dec2hex(value,2); % convert to 2 digit hex
M2_row=hex2dec(hexStr(1))+1;% convert first digit to decimal
M2_col=hex2dec(hexStr(2))+1;% convert first digit to decimal
% %Debugging

% % place the original rows and columns in M2
%
M2(M2_row, M2_col)=row*16+col;
end
end
M2;

If anyone can help.Hello, I have two matrices, M1 ans M2. They both have rows and columns fom 0 to 9 and a to f (16 * 16) matrix.
Mi matrix is filled with random numbers form 0 to 255 in hex decimal notation.

The second matrix M2 is empty but has the same row and column indexs. I need to fill M2 with the help of M1. For eample, in M2, at row 3 and colm 4 ( 34,the first entry of M1) will be 00 ( index of the first entry of M1). In M2 at row 6 col A, the entry will be 01, at row 2 col 9, the entry will be 02, and so on.
%% % Initialize M2 with zeros
M2 = zeros(16,16);
for row = 0:15
for col = 0:15
% get the value from M1
value= M1(row+1,col+1);

% Convert the values to row and col indices for M2

hexStr=dec2hex(value,2); % convert to 2 digit hex
M2_row=hex2dec(hexStr(1))+1;% convert first digit to decimal
M2_col=hex2dec(hexStr(2))+1;% convert first digit to decimal
% %Debugging

% % place the original rows and columns in M2
%
M2(M2_row, M2_col)=row*16+col;
end
end
M2;

If anyone can help. Hello, I have two matrices, M1 ans M2. They both have rows and columns fom 0 to 9 and a to f (16 * 16) matrix.
Mi matrix is filled with random numbers form 0 to 255 in hex decimal notation.

The second matrix M2 is empty but has the same row and column indexs. I need to fill M2 with the help of M1. For eample, in M2, at row 3 and colm 4 ( 34,the first entry of M1) will be 00 ( index of the first entry of M1). In M2 at row 6 col A, the entry will be 01, at row 2 col 9, the entry will be 02, and so on.
%% % Initialize M2 with zeros
M2 = zeros(16,16);
for row = 0:15
for col = 0:15
% get the value from M1
value= M1(row+1,col+1);

% Convert the values to row and col indices for M2

hexStr=dec2hex(value,2); % convert to 2 digit hex
M2_row=hex2dec(hexStr(1))+1;% convert first digit to decimal
M2_col=hex2dec(hexStr(2))+1;% convert first digit to decimal
% %Debugging

% % place the original rows and columns in M2
%
M2(M2_row, M2_col)=row*16+col;
end
end
M2;

If anyone can help. matrix indexes MATLAB Answers — New Questions

​

Hello everyone
Hello everyone,I’m encountering some difficulties while analyzing EEG data using EEGLab. After clustering ICs using the K-means algorithm, I’m unsure about how to select meaningful clusters for further analysis.Specifically, I’d like to plot ERSPs for a chosen cluster (e.g., cluster 2) and then apply bootstrapping to the entire epoch. Subsequently, I want to set all non-significant ERSP values (p > 0.05) compared to this distribution to zero dB.Hello everyone
Hello everyone,I’m encountering some difficulties while analyzing EEG data using EEGLab. After clustering ICs using the K-means algorithm, I’m unsure about how to select meaningful clusters for further analysis.Specifically, I’d like to plot ERSPs for a chosen cluster (e.g., cluster 2) and then apply bootstrapping to the entire epoch. Subsequently, I want to set all non-significant ERSP values (p > 0.05) compared to this distribution to zero dB. Hello everyone
Hello everyone,I’m encountering some difficulties while analyzing EEG data using EEGLab. After clustering ICs using the K-means algorithm, I’m unsure about how to select meaningful clusters for further analysis.Specifically, I’d like to plot ERSPs for a chosen cluster (e.g., cluster 2) and then apply bootstrapping to the entire epoch. Subsequently, I want to set all non-significant ERSP values (p > 0.05) compared to this distribution to zero dB. matlab, eeglab, eeg, ersp, clustering MATLAB Answers — New Questions

​

Hello everyone, I have a question
I have a continuous system that updates every 0.001 seconds.
In addition to this, I have a controller that should be updated with a delay of at least 100 times the system (ie every 100 iterations it is updated)

My question is this:
Do I need to calculate the signal and its derivatives and in addition the error and its derivative with T or with dt? I’m not sure..

The system is denoted as y_pk
The required signal is denoted as y_dk

This is part of my code:
dt=0.001; %discrete time
T=0.1; %continoues time
y_d = sin(1 * time); % desired signal
dy_d = diff(y_d) / dt;
ddy_d = diff(dy_d) / dt;
%initialize:
y_dk1=0; dy_dk1=0; uk2=0; uk1=0;
firstrun=true;
for j = 1:length(time)
y_dk = y_d(j);
if mod(j, 100) == 0 ||firstrun
T=dt; % ??
dy_dk = (y_dk-y_dk1)/T; %derv of y_dk
ddy_dk = (dy_dk – dy_dk1)/T; %second derv of y_dk
ek = y_dk – y_pk;
de_k = dy_dk – dy_pk;

uk = (1/c)*(a1*dy_pk + b1*(dy_pk)^2 + a2*y_pk + b2*(y_pk)^2 +k_p*ek + k_d*de_k + ddy_dk); % u: control input
% Store data for nex step:
uk2=uk1;
uk1=uk;
firstrun=false;
end
y_pk = …. ; %update each 0.001 second (continoues time)
dy_pk = (y_pk-y_pk1)/dt; %derv of y_pk

% Store data for nex step: system
y_pk2=y_pk1;
y_pk1=y_pk;

Thank you all!! :)Hello everyone, I have a question
I have a continuous system that updates every 0.001 seconds.
In addition to this, I have a controller that should be updated with a delay of at least 100 times the system (ie every 100 iterations it is updated)

My question is this:
Do I need to calculate the signal and its derivatives and in addition the error and its derivative with T or with dt? I’m not sure..

The system is denoted as y_pk
The required signal is denoted as y_dk

This is part of my code:
dt=0.001; %discrete time
T=0.1; %continoues time
y_d = sin(1 * time); % desired signal
dy_d = diff(y_d) / dt;
ddy_d = diff(dy_d) / dt;
%initialize:
y_dk1=0; dy_dk1=0; uk2=0; uk1=0;
firstrun=true;
for j = 1:length(time)
y_dk = y_d(j);
if mod(j, 100) == 0 ||firstrun
T=dt; % ??
dy_dk = (y_dk-y_dk1)/T; %derv of y_dk
ddy_dk = (dy_dk – dy_dk1)/T; %second derv of y_dk
ek = y_dk – y_pk;
de_k = dy_dk – dy_pk;

uk = (1/c)*(a1*dy_pk + b1*(dy_pk)^2 + a2*y_pk + b2*(y_pk)^2 +k_p*ek + k_d*de_k + ddy_dk); % u: control input
% Store data for nex step:
uk2=uk1;
uk1=uk;
firstrun=false;
end
y_pk = …. ; %update each 0.001 second (continoues time)
dy_pk = (y_pk-y_pk1)/dt; %derv of y_pk

% Store data for nex step: system
y_pk2=y_pk1;
y_pk1=y_pk;

Thank you all!! 🙂 Hello everyone, I have a question
I have a continuous system that updates every 0.001 seconds.
In addition to this, I have a controller that should be updated with a delay of at least 100 times the system (ie every 100 iterations it is updated)

My question is this:
Do I need to calculate the signal and its derivatives and in addition the error and its derivative with T or with dt? I’m not sure..

The system is denoted as y_pk
The required signal is denoted as y_dk

This is part of my code:
dt=0.001; %discrete time
T=0.1; %continoues time
y_d = sin(1 * time); % desired signal
dy_d = diff(y_d) / dt;
ddy_d = diff(dy_d) / dt;
%initialize:
y_dk1=0; dy_dk1=0; uk2=0; uk1=0;
firstrun=true;
for j = 1:length(time)
y_dk = y_d(j);
if mod(j, 100) == 0 ||firstrun
T=dt; % ??
dy_dk = (y_dk-y_dk1)/T; %derv of y_dk
ddy_dk = (dy_dk – dy_dk1)/T; %second derv of y_dk
ek = y_dk – y_pk;
de_k = dy_dk – dy_pk;

uk = (1/c)*(a1*dy_pk + b1*(dy_pk)^2 + a2*y_pk + b2*(y_pk)^2 +k_p*ek + k_d*de_k + ddy_dk); % u: control input
% Store data for nex step:
uk2=uk1;
uk1=uk;
firstrun=false;
end
y_pk = …. ; %update each 0.001 second (continoues time)
dy_pk = (y_pk-y_pk1)/dt; %derv of y_pk

% Store data for nex step: system
y_pk2=y_pk1;
y_pk1=y_pk;

Thank you all!! 🙂 control, discrete system, controller MATLAB Answers — New Questions

​

Good morning everyone,
I am trying to extract the fracture aperture from a tomography image. I have written the code for the extraction, and everything works using the fitting smoothing spline (createFit function attached) for clu_sli_12.mat and I_12.mat, but not on clu_sli_1400.mat and I_1400.mat in which we have non-unique data (presented directly in the code below). The variable clu_sli_1400 represents the point cloud used as reference for selecting data in the tomography slice (variable I_1400). We tried using other functions in matlab, without success. Can anybody help suggesting a specific function for fitting this data?
thanks in advance for any help!!!!!!!!
clear;
close all;
clc;
warning off;

load(‘I_1400.mat’)
load(‘clu_sli_1400.mat’)

voxel_size=70.69*10^-3;

%% Iteration on the slices
mean_widths = [];
std_widths = [];
min_widths = [];
max_widths = [];
areas = [];

% for ii = 12 %we have to iterate to ii = 12 and ii = 1400 (the last one not working properly on the fitting)

% fitting data
X = clu_sli(:,1);
Y = clu_sli(:,2);

% Create model fitting
[fitresult, gof] = createFit(X, Y);

% Calcolare i valori della curva di fitting
vv = ppval(fitresult.p, X);

% Calcolare la prima derivata della curva di fitting
d1 = ppval(fnder(fitresult.p, 1), X);

% Calcolare la normale alla curva di fitting
epsilon = 1e-10; % Evita la divisione per zero
nn = 1 ./ (d1 + epsilon);

% Calcolare gli angoli delle normali
theta = atan(-1 ./ d1); % Angolo della normale (in radianti)

% Calcolare le componenti della normale
normal_x = cos(theta); % Componente x della normale
normal_y = sin(theta); % Componente y della normale

% Calcolare l’ampiezza della crepa e visualizzare i segmenti
span = 50;
largh = zeros(size(clu_sli, 1), 1);

% Creare una figura con due subplot
figure;

% Primo subplot: Visualizzare la slice con i pallini rossi
subplot(1, 2, 1);
pcolor(I);
shading interp;
axis equal;
hold on;
scatter(clu_sli(:,1), clu_sli(:,2), ‘r’);
xlabel(‘X’);
ylabel(‘Y’);

% Secondo subplot: Visualizzare la slice con i segmenti verdi
subplot(1, 2, 2);
pcolor(I);
shading interp;
axis equal;
hold on;
scatter(clu_sli(:,1), clu_sli(:,2), ‘r’);
plot(X, vv, ‘b’, ‘LineWidth’, 2);
quiver(X, vv, normal_x, normal_y, 0.5, ‘k’, ‘LineWidth’, 1);

intensity_profiles = cell(size(clu_sli, 1), 1);
vet_in = zeros(size(clu_sli, 1), 2);
vet_fin = zeros(size(clu_sli, 1), 2);

for kk = 1:size(clu_sli, 1)
c = round(X(kk));
r = round(vv(kk));

% Genera i pixel lungo la normale
x_profile = round(c + (-span:span) * normal_x(kk));
y_profile = round(r + (-span:span) * normal_y(kk));

% Assicurarsi che i pixel siano all’interno dell’immagine
valid_idx = x_profile > 0 & x_profile <= size(I, 2) & y_profile > 0 & y_profile <= size(I, 1);
x_profile = x_profile(valid_idx);
y_profile = y_profile(valid_idx);

% Estrarre i valori del profilo dall’immagine
profilo = I(sub2ind(size(I), y_profile, x_profile));
intensity_profiles{kk} = profilo;

% Calcolare l’ampiezza della crepa
largh(kk) = conv_crepa(profilo);

% Calcolare i punti del segmento di larghezza
half_width = largh(kk) / 2;
segment_x = [c – half_width * normal_x(kk), c + half_width * normal_x(kk)];
segment_y = [r – half_width * normal_y(kk), r + half_width * normal_y(kk)];

vet_in(kk,1) = c – half_width * normal_x(kk);
vet_fin(kk,1) = c + half_width * normal_x(kk);

vet_in(kk,2) = r – half_width * normal_y(kk);
vet_fin(kk,2) = r + half_width * normal_y(kk);

% Visualizzare il segmento di larghezza
plot(segment_x, segment_y, ‘g’, ‘LineWidth’, 2);
end

hold off;

% Calcolare la media, la deviazione standard, il minimo e il massimo della larghezza della crepa
mean_width = mean(largh, ‘omitnan’);
std_width = std(largh, ‘omitnan’);
min_width = min(largh, [], ‘omitnan’);
max_width = max(largh, [], ‘omitnan’);

mean_widths = [mean_widths; mean_width];
std_widths = [std_widths; std_width];
min_widths = [min_widths; min_width];
max_widths = [max_widths; max_width];

% Calcolo dell’area all’interno delle curve usando trapz
% Creazione del poligono per calcolare l’area
all_x = [vet_in(:,1); flip(vet_fin(:,1))];
all_y = [vet_in(:,2); flip(vet_fin(:,2))];

% Rimuovi i NaN
valid_idx = ~isnan(all_x) & ~isnan(all_y);
all_x = all_x(valid_idx);
all_y = all_y(valid_idx);

% Assicurarsi che i punti formino un poligono chiuso
if length(all_x) > 1 && (all_x(1) ~= all_x(end) || all_y(1) ~= all_y(end))
all_x(end+1) = all_x(1);
all_y(end+1) = all_y(1);
end

% Calcolo dell’area del poligono usando trapz
area = trapz(all_x, all_y);
areas = [areas; area];

title(‘Area Between Curves per Slice’);

%% calcoli in mm

areas_mm = areas.* voxel_size^2
std_mm=std_widths.*voxel_size
mean_mm=mean_widths.*voxel_size
min_mm=min_widths.*voxel_size
max_mm=max_widths.*voxel_sizeGood morning everyone,
I am trying to extract the fracture aperture from a tomography image. I have written the code for the extraction, and everything works using the fitting smoothing spline (createFit function attached) for clu_sli_12.mat and I_12.mat, but not on clu_sli_1400.mat and I_1400.mat in which we have non-unique data (presented directly in the code below). The variable clu_sli_1400 represents the point cloud used as reference for selecting data in the tomography slice (variable I_1400). We tried using other functions in matlab, without success. Can anybody help suggesting a specific function for fitting this data?
thanks in advance for any help!!!!!!!!
clear;
close all;
clc;
warning off;

load(‘I_1400.mat’)
load(‘clu_sli_1400.mat’)

voxel_size=70.69*10^-3;

%% Iteration on the slices
mean_widths = [];
std_widths = [];
min_widths = [];
max_widths = [];
areas = [];

% for ii = 12 %we have to iterate to ii = 12 and ii = 1400 (the last one not working properly on the fitting)

% fitting data
X = clu_sli(:,1);
Y = clu_sli(:,2);

% Create model fitting
[fitresult, gof] = createFit(X, Y);

% Calcolare i valori della curva di fitting
vv = ppval(fitresult.p, X);

% Calcolare la prima derivata della curva di fitting
d1 = ppval(fnder(fitresult.p, 1), X);

% Calcolare la normale alla curva di fitting
epsilon = 1e-10; % Evita la divisione per zero
nn = 1 ./ (d1 + epsilon);

% Calcolare gli angoli delle normali
theta = atan(-1 ./ d1); % Angolo della normale (in radianti)

% Calcolare le componenti della normale
normal_x = cos(theta); % Componente x della normale
normal_y = sin(theta); % Componente y della normale

% Calcolare l’ampiezza della crepa e visualizzare i segmenti
span = 50;
largh = zeros(size(clu_sli, 1), 1);

% Creare una figura con due subplot
figure;

% Primo subplot: Visualizzare la slice con i pallini rossi
subplot(1, 2, 1);
pcolor(I);
shading interp;
axis equal;
hold on;
scatter(clu_sli(:,1), clu_sli(:,2), ‘r’);
xlabel(‘X’);
ylabel(‘Y’);

% Secondo subplot: Visualizzare la slice con i segmenti verdi
subplot(1, 2, 2);
pcolor(I);
shading interp;
axis equal;
hold on;
scatter(clu_sli(:,1), clu_sli(:,2), ‘r’);
plot(X, vv, ‘b’, ‘LineWidth’, 2);
quiver(X, vv, normal_x, normal_y, 0.5, ‘k’, ‘LineWidth’, 1);

intensity_profiles = cell(size(clu_sli, 1), 1);
vet_in = zeros(size(clu_sli, 1), 2);
vet_fin = zeros(size(clu_sli, 1), 2);

for kk = 1:size(clu_sli, 1)
c = round(X(kk));
r = round(vv(kk));

% Genera i pixel lungo la normale
x_profile = round(c + (-span:span) * normal_x(kk));
y_profile = round(r + (-span:span) * normal_y(kk));

% Assicurarsi che i pixel siano all’interno dell’immagine
valid_idx = x_profile > 0 & x_profile <= size(I, 2) & y_profile > 0 & y_profile <= size(I, 1);
x_profile = x_profile(valid_idx);
y_profile = y_profile(valid_idx);

% Estrarre i valori del profilo dall’immagine
profilo = I(sub2ind(size(I), y_profile, x_profile));
intensity_profiles{kk} = profilo;

% Calcolare l’ampiezza della crepa
largh(kk) = conv_crepa(profilo);

% Calcolare i punti del segmento di larghezza
half_width = largh(kk) / 2;
segment_x = [c – half_width * normal_x(kk), c + half_width * normal_x(kk)];
segment_y = [r – half_width * normal_y(kk), r + half_width * normal_y(kk)];

vet_in(kk,1) = c – half_width * normal_x(kk);
vet_fin(kk,1) = c + half_width * normal_x(kk);

vet_in(kk,2) = r – half_width * normal_y(kk);
vet_fin(kk,2) = r + half_width * normal_y(kk);

% Visualizzare il segmento di larghezza
plot(segment_x, segment_y, ‘g’, ‘LineWidth’, 2);
end

hold off;

% Calcolare la media, la deviazione standard, il minimo e il massimo della larghezza della crepa
mean_width = mean(largh, ‘omitnan’);
std_width = std(largh, ‘omitnan’);
min_width = min(largh, [], ‘omitnan’);
max_width = max(largh, [], ‘omitnan’);

mean_widths = [mean_widths; mean_width];
std_widths = [std_widths; std_width];
min_widths = [min_widths; min_width];
max_widths = [max_widths; max_width];

% Calcolo dell’area all’interno delle curve usando trapz
% Creazione del poligono per calcolare l’area
all_x = [vet_in(:,1); flip(vet_fin(:,1))];
all_y = [vet_in(:,2); flip(vet_fin(:,2))];

% Rimuovi i NaN
valid_idx = ~isnan(all_x) & ~isnan(all_y);
all_x = all_x(valid_idx);
all_y = all_y(valid_idx);

% Assicurarsi che i punti formino un poligono chiuso
if length(all_x) > 1 && (all_x(1) ~= all_x(end) || all_y(1) ~= all_y(end))
all_x(end+1) = all_x(1);
all_y(end+1) = all_y(1);
end

% Calcolo dell’area del poligono usando trapz
area = trapz(all_x, all_y);
areas = [areas; area];

title(‘Area Between Curves per Slice’);

%% calcoli in mm

areas_mm = areas.* voxel_size^2
std_mm=std_widths.*voxel_size
mean_mm=mean_widths.*voxel_size
min_mm=min_widths.*voxel_size
max_mm=max_widths.*voxel_size Good morning everyone,
I am trying to extract the fracture aperture from a tomography image. I have written the code for the extraction, and everything works using the fitting smoothing spline (createFit function attached) for clu_sli_12.mat and I_12.mat, but not on clu_sli_1400.mat and I_1400.mat in which we have non-unique data (presented directly in the code below). The variable clu_sli_1400 represents the point cloud used as reference for selecting data in the tomography slice (variable I_1400). We tried using other functions in matlab, without success. Can anybody help suggesting a specific function for fitting this data?
thanks in advance for any help!!!!!!!!
clear;
close all;
clc;
warning off;

load(‘I_1400.mat’)
load(‘clu_sli_1400.mat’)

voxel_size=70.69*10^-3;

%% Iteration on the slices
mean_widths = [];
std_widths = [];
min_widths = [];
max_widths = [];
areas = [];

% for ii = 12 %we have to iterate to ii = 12 and ii = 1400 (the last one not working properly on the fitting)

% fitting data
X = clu_sli(:,1);
Y = clu_sli(:,2);

% Create model fitting
[fitresult, gof] = createFit(X, Y);

% Calcolare i valori della curva di fitting
vv = ppval(fitresult.p, X);

% Calcolare la prima derivata della curva di fitting
d1 = ppval(fnder(fitresult.p, 1), X);

% Calcolare la normale alla curva di fitting
epsilon = 1e-10; % Evita la divisione per zero
nn = 1 ./ (d1 + epsilon);

% Calcolare gli angoli delle normali
theta = atan(-1 ./ d1); % Angolo della normale (in radianti)

% Calcolare le componenti della normale
normal_x = cos(theta); % Componente x della normale
normal_y = sin(theta); % Componente y della normale

% Calcolare l’ampiezza della crepa e visualizzare i segmenti
span = 50;
largh = zeros(size(clu_sli, 1), 1);

% Creare una figura con due subplot
figure;

% Primo subplot: Visualizzare la slice con i pallini rossi
subplot(1, 2, 1);
pcolor(I);
shading interp;
axis equal;
hold on;
scatter(clu_sli(:,1), clu_sli(:,2), ‘r’);
xlabel(‘X’);
ylabel(‘Y’);

% Secondo subplot: Visualizzare la slice con i segmenti verdi
subplot(1, 2, 2);
pcolor(I);
shading interp;
axis equal;
hold on;
scatter(clu_sli(:,1), clu_sli(:,2), ‘r’);
plot(X, vv, ‘b’, ‘LineWidth’, 2);
quiver(X, vv, normal_x, normal_y, 0.5, ‘k’, ‘LineWidth’, 1);

intensity_profiles = cell(size(clu_sli, 1), 1);
vet_in = zeros(size(clu_sli, 1), 2);
vet_fin = zeros(size(clu_sli, 1), 2);

for kk = 1:size(clu_sli, 1)
c = round(X(kk));
r = round(vv(kk));

% Genera i pixel lungo la normale
x_profile = round(c + (-span:span) * normal_x(kk));
y_profile = round(r + (-span:span) * normal_y(kk));

% Assicurarsi che i pixel siano all’interno dell’immagine
valid_idx = x_profile > 0 & x_profile <= size(I, 2) & y_profile > 0 & y_profile <= size(I, 1);
x_profile = x_profile(valid_idx);
y_profile = y_profile(valid_idx);

% Estrarre i valori del profilo dall’immagine
profilo = I(sub2ind(size(I), y_profile, x_profile));
intensity_profiles{kk} = profilo;

% Calcolare l’ampiezza della crepa
largh(kk) = conv_crepa(profilo);

% Calcolare i punti del segmento di larghezza
half_width = largh(kk) / 2;
segment_x = [c – half_width * normal_x(kk), c + half_width * normal_x(kk)];
segment_y = [r – half_width * normal_y(kk), r + half_width * normal_y(kk)];

vet_in(kk,1) = c – half_width * normal_x(kk);
vet_fin(kk,1) = c + half_width * normal_x(kk);

vet_in(kk,2) = r – half_width * normal_y(kk);
vet_fin(kk,2) = r + half_width * normal_y(kk);

% Visualizzare il segmento di larghezza
plot(segment_x, segment_y, ‘g’, ‘LineWidth’, 2);
end

hold off;

% Calcolare la media, la deviazione standard, il minimo e il massimo della larghezza della crepa
mean_width = mean(largh, ‘omitnan’);
std_width = std(largh, ‘omitnan’);
min_width = min(largh, [], ‘omitnan’);
max_width = max(largh, [], ‘omitnan’);

mean_widths = [mean_widths; mean_width];
std_widths = [std_widths; std_width];
min_widths = [min_widths; min_width];
max_widths = [max_widths; max_width];

% Calcolo dell’area all’interno delle curve usando trapz
% Creazione del poligono per calcolare l’area
all_x = [vet_in(:,1); flip(vet_fin(:,1))];
all_y = [vet_in(:,2); flip(vet_fin(:,2))];

% Rimuovi i NaN
valid_idx = ~isnan(all_x) & ~isnan(all_y);
all_x = all_x(valid_idx);
all_y = all_y(valid_idx);

% Assicurarsi che i punti formino un poligono chiuso
if length(all_x) > 1 && (all_x(1) ~= all_x(end) || all_y(1) ~= all_y(end))
all_x(end+1) = all_x(1);
all_y(end+1) = all_y(1);
end

% Calcolo dell’area del poligono usando trapz
area = trapz(all_x, all_y);
areas = [areas; area];

title(‘Area Between Curves per Slice’);

%% calcoli in mm

areas_mm = areas.* voxel_size^2
std_mm=std_widths.*voxel_size
mean_mm=mean_widths.*voxel_size
min_mm=min_widths.*voxel_size
max_mm=max_widths.*voxel_size curve fitting, non-unique data, createfit MATLAB Answers — New Questions

​

I have a piece of code that is currently calculating stress/strain and force/extension.
It groups them into catagories based on the temperature and time they were sintered for.
I would like to add histograsm that can also show the modulus of elasticity and the ultimate tensile strength for these groups. The modulus of elasticity should be for say a range from ~0.05 to 0.15 extension, or similar. The error bars should show +/- the stdev for the data sets
would also like to make it so that the legend expressly states the sintering cycle.
For example ‘conventional sintering’ two-stage sintering’, ‘sintered at x°C for x hours’ etc
below is the code
clear
close all

Nsamples = 14;
G=ones(1,2*Nsamples);
G(1:12)=1;
G([2:1:4 10:1:11])=2;
G(5:1:9)=3;
G(13:14)=4;
%%
load(‘E:DICsintering_dataOXmcdata.mat’)

% G(15)=5;
% G(16)=6;
% G(13:14)=5;
% M={‘x’,’+’,’*’,’s’,’d’}; % a marker for each group
C={‘r’,’b’,’m’,’g’,’m’,’y’}; % and a color
L = 90; % coupon gauge section (mm)
A = 18; % coupon cross-sectional area (mm^2)
figure(1);
hold on
figure(2);
hold on
for k = 1:Nsamples
file_name = [‘Sample’ num2str(k) ‘.csv’];
T = readtable(file_name,’VariableNamingRule’,’preserve’);
[test(k).disp,test(k).force,test(k).eyy] = post_fun(T);
ig=G(k); % get group id for kth file
figure(1)
% h(k) = plot(test(k).disp,test(k).force,’LineWidth’,2,’Marker’,M(ig),’Color’,C{ig});
h(k) = plot(test(k).disp,test(k).force,’LineWidth’,2,’Color’,C{ig});
figure(2)
% g(k) = plot(test(k).eyy*100,1000*test(k).force/A,’LineWidth’,2,’Marker’,M(ig),’Color’,C{ig});
g(k) = plot(test(k).eyy*100,1000*test(k).force/A,’LineWidth’,2,’Color’,C{ig});
leg_string{k} = [‘Sintering schedule ‘ num2str(ig)];
end
%
%%
figure(1)
set(gcf,’Position’,[200 200 1024 768],’Color’,[1 1 1]);
set(gca,’FontSize’,24);
xlabel(‘Displacement (mm)’,’FontSize’,32,’Interpreter’,’latex’);
ylabel(‘Force (kN)’,’FontSize’,32,’Interpreter’,’latex’);
box on
grid on
title(‘Force / extension curves for of various sintering schedules’,’FontSize’, 20,’Interpreter’,’latex’)
legend(h([1 2 5 13]),leg_string([1 2 5 13]),’FontSize’,28,’Interpreter’,’latex’,’location’,’southeast’);
%%
figure(2)
set(gcf,’Position’,[200 200 1024 768],’Color’,[1 1 1]);
set(gca,’FontSize’,24);
xlabel(‘Strain (%)’,’FontSize’,32,’Interpreter’,’latex’);
ylabel(‘Stress (MPa)’,’FontSize’,32,’Interpreter’,’latex’);
box on
grid on
title(‘Stress / strain curves for of various sintering schedules’,’FontSize’, 20,’Interpreter’,’latex’)
legend(g([1 2 5 13]),leg_string([1 2 5 13]),’FontSize’,28,’Interpreter’,’latex’,’location’,’southeast’);
%
plot(OXmcdata(:,1),OXmcdata(:,2),’DisplayName’,’Farhandi et al (2021)’,’LineWidth’,3,’Color’,’k’,’LineStyle’,’–‘);
%
xlim([0 3E-1]);
ylim([0 350]);
%%

function [displ,force,eyy] = post_fun(T)
displ=T{:,4};
force=T{:,end};
eyy =T{:,10};
ix=all(isfinite([displ force eyy]),2) & all([displ force eyy]>0,2);
displ=displ(ix);
force=force(ix);
eyy=eyy(ix);
[~,imax] = max(force);
displ(imax+1:end) = [];
force(imax+1:end) = [];
eyy(imax+1:end) = [];
end
Hopefully you can help.

Thanks

AlexI have a piece of code that is currently calculating stress/strain and force/extension.
It groups them into catagories based on the temperature and time they were sintered for.
I would like to add histograsm that can also show the modulus of elasticity and the ultimate tensile strength for these groups. The modulus of elasticity should be for say a range from ~0.05 to 0.15 extension, or similar. The error bars should show +/- the stdev for the data sets
would also like to make it so that the legend expressly states the sintering cycle.
For example ‘conventional sintering’ two-stage sintering’, ‘sintered at x°C for x hours’ etc
below is the code
clear
close all

Nsamples = 14;
G=ones(1,2*Nsamples);
G(1:12)=1;
G([2:1:4 10:1:11])=2;
G(5:1:9)=3;
G(13:14)=4;
%%
load(‘E:DICsintering_dataOXmcdata.mat’)

% G(15)=5;
% G(16)=6;
% G(13:14)=5;
% M={‘x’,’+’,’*’,’s’,’d’}; % a marker for each group
C={‘r’,’b’,’m’,’g’,’m’,’y’}; % and a color
L = 90; % coupon gauge section (mm)
A = 18; % coupon cross-sectional area (mm^2)
figure(1);
hold on
figure(2);
hold on
for k = 1:Nsamples
file_name = [‘Sample’ num2str(k) ‘.csv’];
T = readtable(file_name,’VariableNamingRule’,’preserve’);
[test(k).disp,test(k).force,test(k).eyy] = post_fun(T);
ig=G(k); % get group id for kth file
figure(1)
% h(k) = plot(test(k).disp,test(k).force,’LineWidth’,2,’Marker’,M(ig),’Color’,C{ig});
h(k) = plot(test(k).disp,test(k).force,’LineWidth’,2,’Color’,C{ig});
figure(2)
% g(k) = plot(test(k).eyy*100,1000*test(k).force/A,’LineWidth’,2,’Marker’,M(ig),’Color’,C{ig});
g(k) = plot(test(k).eyy*100,1000*test(k).force/A,’LineWidth’,2,’Color’,C{ig});
leg_string{k} = [‘Sintering schedule ‘ num2str(ig)];
end
%
%%
figure(1)
set(gcf,’Position’,[200 200 1024 768],’Color’,[1 1 1]);
set(gca,’FontSize’,24);
xlabel(‘Displacement (mm)’,’FontSize’,32,’Interpreter’,’latex’);
ylabel(‘Force (kN)’,’FontSize’,32,’Interpreter’,’latex’);
box on
grid on
title(‘Force / extension curves for of various sintering schedules’,’FontSize’, 20,’Interpreter’,’latex’)
legend(h([1 2 5 13]),leg_string([1 2 5 13]),’FontSize’,28,’Interpreter’,’latex’,’location’,’southeast’);
%%
figure(2)
set(gcf,’Position’,[200 200 1024 768],’Color’,[1 1 1]);
set(gca,’FontSize’,24);
xlabel(‘Strain (%)’,’FontSize’,32,’Interpreter’,’latex’);
ylabel(‘Stress (MPa)’,’FontSize’,32,’Interpreter’,’latex’);
box on
grid on
title(‘Stress / strain curves for of various sintering schedules’,’FontSize’, 20,’Interpreter’,’latex’)
legend(g([1 2 5 13]),leg_string([1 2 5 13]),’FontSize’,28,’Interpreter’,’latex’,’location’,’southeast’);
%
plot(OXmcdata(:,1),OXmcdata(:,2),’DisplayName’,’Farhandi et al (2021)’,’LineWidth’,3,’Color’,’k’,’LineStyle’,’–‘);
%
xlim([0 3E-1]);
ylim([0 350]);
%%

function [displ,force,eyy] = post_fun(T)
displ=T{:,4};
force=T{:,end};
eyy =T{:,10};
ix=all(isfinite([displ force eyy]),2) & all([displ force eyy]>0,2);
displ=displ(ix);
force=force(ix);
eyy=eyy(ix);
[~,imax] = max(force);
displ(imax+1:end) = [];
force(imax+1:end) = [];
eyy(imax+1:end) = [];
end
Hopefully you can help.

Thanks

Alex I have a piece of code that is currently calculating stress/strain and force/extension.
It groups them into catagories based on the temperature and time they were sintered for.
I would like to add histograsm that can also show the modulus of elasticity and the ultimate tensile strength for these groups. The modulus of elasticity should be for say a range from ~0.05 to 0.15 extension, or similar. The error bars should show +/- the stdev for the data sets
would also like to make it so that the legend expressly states the sintering cycle.
For example ‘conventional sintering’ two-stage sintering’, ‘sintered at x°C for x hours’ etc
below is the code
clear
close all

Nsamples = 14;
G=ones(1,2*Nsamples);
G(1:12)=1;
G([2:1:4 10:1:11])=2;
G(5:1:9)=3;
G(13:14)=4;
%%
load(‘E:DICsintering_dataOXmcdata.mat’)

% G(15)=5;
% G(16)=6;
% G(13:14)=5;
% M={‘x’,’+’,’*’,’s’,’d’}; % a marker for each group
C={‘r’,’b’,’m’,’g’,’m’,’y’}; % and a color
L = 90; % coupon gauge section (mm)
A = 18; % coupon cross-sectional area (mm^2)
figure(1);
hold on
figure(2);
hold on
for k = 1:Nsamples
file_name = [‘Sample’ num2str(k) ‘.csv’];
T = readtable(file_name,’VariableNamingRule’,’preserve’);
[test(k).disp,test(k).force,test(k).eyy] = post_fun(T);
ig=G(k); % get group id for kth file
figure(1)
% h(k) = plot(test(k).disp,test(k).force,’LineWidth’,2,’Marker’,M(ig),’Color’,C{ig});
h(k) = plot(test(k).disp,test(k).force,’LineWidth’,2,’Color’,C{ig});
figure(2)
% g(k) = plot(test(k).eyy*100,1000*test(k).force/A,’LineWidth’,2,’Marker’,M(ig),’Color’,C{ig});
g(k) = plot(test(k).eyy*100,1000*test(k).force/A,’LineWidth’,2,’Color’,C{ig});
leg_string{k} = [‘Sintering schedule ‘ num2str(ig)];
end
%
%%
figure(1)
set(gcf,’Position’,[200 200 1024 768],’Color’,[1 1 1]);
set(gca,’FontSize’,24);
xlabel(‘Displacement (mm)’,’FontSize’,32,’Interpreter’,’latex’);
ylabel(‘Force (kN)’,’FontSize’,32,’Interpreter’,’latex’);
box on
grid on
title(‘Force / extension curves for of various sintering schedules’,’FontSize’, 20,’Interpreter’,’latex’)
legend(h([1 2 5 13]),leg_string([1 2 5 13]),’FontSize’,28,’Interpreter’,’latex’,’location’,’southeast’);
%%
figure(2)
set(gcf,’Position’,[200 200 1024 768],’Color’,[1 1 1]);
set(gca,’FontSize’,24);
xlabel(‘Strain (%)’,’FontSize’,32,’Interpreter’,’latex’);
ylabel(‘Stress (MPa)’,’FontSize’,32,’Interpreter’,’latex’);
box on
grid on
title(‘Stress / strain curves for of various sintering schedules’,’FontSize’, 20,’Interpreter’,’latex’)
legend(g([1 2 5 13]),leg_string([1 2 5 13]),’FontSize’,28,’Interpreter’,’latex’,’location’,’southeast’);
%
plot(OXmcdata(:,1),OXmcdata(:,2),’DisplayName’,’Farhandi et al (2021)’,’LineWidth’,3,’Color’,’k’,’LineStyle’,’–‘);
%
xlim([0 3E-1]);
ylim([0 350]);
%%

function [displ,force,eyy] = post_fun(T)
displ=T{:,4};
force=T{:,end};
eyy =T{:,10};
ix=all(isfinite([displ force eyy]),2) & all([displ force eyy]>0,2);
displ=displ(ix);
force=force(ix);
eyy=eyy(ix);
[~,imax] = max(force);
displ(imax+1:end) = [];
force(imax+1:end) = [];
eyy(imax+1:end) = [];
end
Hopefully you can help.

Thanks

Alex histogram, bar, error, errorbars MATLAB Answers — New Questions

​

My function is . I used following code to interpolation. How can I find the coefficients of the interpolated function ?
Eta = load(‘Eta.txt’);
t = 0:5:(8.25*60); % time
x = [0,3750,7500,8000]; % distance
[X,T] = meshgrid(x,t);

% interpolation
[xq,tq] = meshgrid(0:2:8000,0:(8.25*60)) ;
Eta_intp = interp2(X,T,eta,xq,tq,’spline’);My function is . I used following code to interpolation. How can I find the coefficients of the interpolated function ?
Eta = load(‘Eta.txt’);
t = 0:5:(8.25*60); % time
x = [0,3750,7500,8000]; % distance
[X,T] = meshgrid(x,t);

% interpolation
[xq,tq] = meshgrid(0:2:8000,0:(8.25*60)) ;
Eta_intp = interp2(X,T,eta,xq,tq,’spline’); My function is . I used following code to interpolation. How can I find the coefficients of the interpolated function ?
Eta = load(‘Eta.txt’);
t = 0:5:(8.25*60); % time
x = [0,3750,7500,8000]; % distance
[X,T] = meshgrid(x,t);

% interpolation
[xq,tq] = meshgrid(0:2:8000,0:(8.25*60)) ;
Eta_intp = interp2(X,T,eta,xq,tq,’spline’); 2d interpolation MATLAB Answers — New Questions

​

Hello,
I imported 2 URDF files using the importrobot function. Then, I created a robotScenario by adding a robotPlatform for each robot. I want to check during the scenario that there are no collisions between the robots.
I tried using the checkCollision function between the 2 robotPlatforms, but this function does not accept non-rigidBodyTree-based platforms.
Then, I tried to create a robotPlatform ‘hand’ from an STL file and attach ‘hand’ to the ‘robot2Scn’ platform using the attach function, but the position of ‘hand’ does not update when I use the checkCollision(robotScn, ‘hand’, IgnoreSelfCollision="on") function. This function returns 1 when ‘robotScn’ is at the initial position of ‘hand’, otherwise it returns 0.
Is there a function or have you a suggestion to check for collisions between 2 robotPlatforms in a robotScenario?
Thank you.
robot = importrobot("robot.urdf");
robot2 = importrobot("robot_2.urdf");

scenario = robotScenario(UpdateRate=1,StopTime=10);
robotScn = robotPlatform("robot1",scenario,RigidBodyTree=robot);
robotScn2 = robotPlatform("robot2",scenario,RigidBodyTree=robot2);

% Trying to manually add a collision zone (Test with and without using clearCollision on robot2)
stlData = stlread(‘hand_link.stl’);
vertices = stlData.Points;
faces = stlData.ConnectivityList;
handMesh = collisionMesh(vertices);
base = robotPlatform("hand",scenario,Collision=handMesh, initialBasePosition[0.2 0.2 0.2]);
attach(robotScn2,"hand","hand_link",ChildToParentTransform=trvec2tform([0 0 0]))
checkCollision(robotScn,"hand", IgnoreSelfCollision="on")Hello,
I imported 2 URDF files using the importrobot function. Then, I created a robotScenario by adding a robotPlatform for each robot. I want to check during the scenario that there are no collisions between the robots.
I tried using the checkCollision function between the 2 robotPlatforms, but this function does not accept non-rigidBodyTree-based platforms.
Then, I tried to create a robotPlatform ‘hand’ from an STL file and attach ‘hand’ to the ‘robot2Scn’ platform using the attach function, but the position of ‘hand’ does not update when I use the checkCollision(robotScn, ‘hand’, IgnoreSelfCollision="on") function. This function returns 1 when ‘robotScn’ is at the initial position of ‘hand’, otherwise it returns 0.
Is there a function or have you a suggestion to check for collisions between 2 robotPlatforms in a robotScenario?
Thank you.
robot = importrobot("robot.urdf");
robot2 = importrobot("robot_2.urdf");

scenario = robotScenario(UpdateRate=1,StopTime=10);
robotScn = robotPlatform("robot1",scenario,RigidBodyTree=robot);
robotScn2 = robotPlatform("robot2",scenario,RigidBodyTree=robot2);

% Trying to manually add a collision zone (Test with and without using clearCollision on robot2)
stlData = stlread(‘hand_link.stl’);
vertices = stlData.Points;
faces = stlData.ConnectivityList;
handMesh = collisionMesh(vertices);
base = robotPlatform("hand",scenario,Collision=handMesh, initialBasePosition[0.2 0.2 0.2]);
attach(robotScn2,"hand","hand_link",ChildToParentTransform=trvec2tform([0 0 0]))
checkCollision(robotScn,"hand", IgnoreSelfCollision="on") Hello,
I imported 2 URDF files using the importrobot function. Then, I created a robotScenario by adding a robotPlatform for each robot. I want to check during the scenario that there are no collisions between the robots.
I tried using the checkCollision function between the 2 robotPlatforms, but this function does not accept non-rigidBodyTree-based platforms.
Then, I tried to create a robotPlatform ‘hand’ from an STL file and attach ‘hand’ to the ‘robot2Scn’ platform using the attach function, but the position of ‘hand’ does not update when I use the checkCollision(robotScn, ‘hand’, IgnoreSelfCollision="on") function. This function returns 1 when ‘robotScn’ is at the initial position of ‘hand’, otherwise it returns 0.
Is there a function or have you a suggestion to check for collisions between 2 robotPlatforms in a robotScenario?
Thank you.
robot = importrobot("robot.urdf");
robot2 = importrobot("robot_2.urdf");

scenario = robotScenario(UpdateRate=1,StopTime=10);
robotScn = robotPlatform("robot1",scenario,RigidBodyTree=robot);
robotScn2 = robotPlatform("robot2",scenario,RigidBodyTree=robot2);

% Trying to manually add a collision zone (Test with and without using clearCollision on robot2)
stlData = stlread(‘hand_link.stl’);
vertices = stlData.Points;
faces = stlData.ConnectivityList;
handMesh = collisionMesh(vertices);
base = robotPlatform("hand",scenario,Collision=handMesh, initialBasePosition[0.2 0.2 0.2]);
attach(robotScn2,"hand","hand_link",ChildToParentTransform=trvec2tform([0 0 0]))
checkCollision(robotScn,"hand", IgnoreSelfCollision="on") robotics, collision MATLAB Answers — New Questions

​

I have a mat file which has a row vector u and a matrix temp_gbo. The size of this matrix is 100 x 4. I want the following:
Each row of the matrix temp_gbo must have the same arrangment of elements as is in u.The mat file and the script "swapping" are attached but the script doesn’t fulfill my requirement. How should I modify this code?I have a mat file which has a row vector u and a matrix temp_gbo. The size of this matrix is 100 x 4. I want the following:
Each row of the matrix temp_gbo must have the same arrangment of elements as is in u.The mat file and the script "swapping" are attached but the script doesn’t fulfill my requirement. How should I modify this code? I have a mat file which has a row vector u and a matrix temp_gbo. The size of this matrix is 100 x 4. I want the following:
Each row of the matrix temp_gbo must have the same arrangment of elements as is in u.The mat file and the script "swapping" are attached but the script doesn’t fulfill my requirement. How should I modify this code? arrangement of elements, row vector, matrix, swapping MATLAB Answers — New Questions

​

I’m a student, I have to build PCA from scratch using Matlab on iris data.
Iris data have 4 features i want to reduce them to 2.
I reached the sorting of eigenvalues step. What is the next step?I’m a student, I have to build PCA from scratch using Matlab on iris data.
Iris data have 4 features i want to reduce them to 2.
I reached the sorting of eigenvalues step. What is the next step? I’m a student, I have to build PCA from scratch using Matlab on iris data.
Iris data have 4 features i want to reduce them to 2.
I reached the sorting of eigenvalues step. What is the next step? pca, features reduction MATLAB Answers — New Questions

​