Tag Archives: matlab
Incorrect answer using ode15s (and others) solving ode with time varying coefficients (rate equations); issue with interpolation?
Hello, I am not getting the correct answer when computing the ode for this set of rate equations. I believe the issue comes with the interpolation step because when I replace the coefficients with the non-interpolated values (A_col, B_col, C_col, D_col in place of aa, bb, cc, dd), I get the correct answer (a lorentzian). I am only using an array of ones as my "time dependent" part of the equation as I am using this as a test to see if I can get the correct answer before I make things more complicated.
I have tried all of the different interpolation schemes available (linear’, ‘nearest’, etc.). I have tried lowering the tolerance.
Any ideas?!
Thanks in advance!!
%clear;
%% ——————- CONSTANTS ——————- %%
tic
m_ArII = 6.6335e-26; %[kg]
e = 1.602e-19; %[C]
c = 299792458; %[m/s]
k_B = 1.380649e-23; %[J/K]
lambda_0 = 6.116616e-7;%[m]
freq_0 = c/lambda_0; %[Hz]
I_0 = 0.5e+9; %[W/m^2]
lineWidth = 1e+9; % [Hz]
A_1_0 = 0.212e+8; %[Hz]
A_1_2 = 0.759e+8; %[Hz]
B_1_0 = 121.9e+11; %[m^2/(J s)]
B_0_1 = 97.55e+11; %[m^2/(J s)]
solidAngle = sin(15*pi/180);
pulseWidth = 6e-9; % [s]
tau_0 = 1e-3;
tau_1 = 8.51e-9;
n_0 = 5e+13; %[1/m^3]
T_i = 0.1; %[eV]
p = 20; %number of intergrals
% —————————————————— %
%———- Limits for integration and diffeqs———-%
freqRange = 30e+9;
freqLimits = linspace(freq_0 – (freqRange*3),freq_0 + (freqRange*3),p); %[Hz]
centralLaserFreq = linspace((freq_0 – freqRange), (freq_0 + freqRange) ,p); %[Hz]
vLimits = linspace(-1e+4,1e+4,p)’; %[m/s]
tspan = [0 pulseWidth];
time_interp = linspace(0,pulseWidth,p)’;
%———————————————%
detuning = (centralLaserFreq – freq_0)./(1e+9);
% preallocated arrays
N_obs_2 = zeros(p,1);
f1intStore = zeros(p,p);
f_1_compare = zeros(p,p*p);
protoStore = zeros(p,p);
%——set up mesgrid for integration——-%
%rows %cols
[F , V] = meshgrid(freqLimits,vLimits);
%——————————————-%
%ion absorption spectrum
L = (A_1_0/pi)./((F – freq_0*(1 + (V/c)) ).^2 + (A_1_0)^2); %freqLimits is a (p x p) array
%———— INITIAL CONDITIONS ————%
f_0_initial = n_0*sqrt(m_ArII/(2*pi*e*T_i)).*exp(-(m_ArII*(vLimits).^2)./(2*e*T_i)); % (p x 1)
beta = (1e-7);
n_1 = beta*n_0;
f_1_initial = n_1*sqrt(m_ArII/(2*pi*e*T_i)).*exp(-(m_ArII*(vLimits).^2)./(2*e*T_i)); % (p x 1);
initialConds = [f_0_initial; f_1_initial]; %(2*p x 1)
%——————————————–%
%time component
theta = ones(p,1);% this is just an array of 1s for now, but will be a temporally dependent term in the future
%loop through different values of centralLaserFreq
for k = 1:p
%laser intensity spectrum
I_of_nu = I_0.*((lineWidth/(2*pi)))./((F – (centralLaserFreq(k) )).^2 + (lineWidth/2)^2);%(p x p)
%effective laser intensity
protoPhi = (1/(4*pi)).*I_of_nu.*L; %(p x p)
%integrate of frequency
protoPhiInt = trapz(freqLimits,protoPhi,2); %(p x 1)
protoStore(:,k) = protoPhiInt;
% (1 x p) (p x 1)
Phi_of_v_t = theta’.*protoPhiInt; %(p x p) % <— Now that the integration over nu is complete, we can bring in the theta term. Didn’t add it previously because it would’ve made things more complicated.
%example of Phi_of_v_t:
% [t1*phi1, t2*phi1, t3*phi1]
% [t1*phi2, t2*phi2, t3*phi2]
% [t1*phi3, t2*phi3, t3*phi3]
%Coeffitients to the rate equations
A = ((1/tau_0) + B_0_1.*Phi_of_v_t); %(p x p)
B = A_1_0 + B_1_0.*Phi_of_v_t; %(p x p)
C = B_0_1.*Phi_of_v_t; %(p x p)
D = ((1/tau_1) + B_1_0.*Phi_of_v_t); %(p x p)
%loop through columns of the coefficients to make things faster than going
%through each element individually
f1intStore = zeros(p,p);
for col = 1:p
[f_1_Sol,tspanOut,f_Sol] = solveLIFRateEqns(A(:,col),B(:,col),C(:,col),D(:,col),initialConds,tspan,time_interp,p);
f_1_timeInt = trapz(tspanOut,f_1_Sol,1); %integrate over time (1 x p)
f1intStore(:,col) = f_1_timeInt; % each column of f1solStore is the f_1 for a different time
end
N_obs_int = (solidAngle/(4*pi)).*A_1_2.*f1intStore(:,1); %(p x 1)
fxintAll2 = trapz(vLimits,N_obs_int); %intgrate over veocity
N_obs_2(k,1) = fxintAll2; %store value
disp(k)
end
%% PLOTTING
[GaussFit, GaussGof] = Gauss2(detuning’, N_obs_2);
% MATLAB’s "gauss1" fitting defines the Gaussian as a1*exp(-((x-b1)/c1)^2) and gives c1 as a result of fitting.
% I’m using a Gaussian defined as a*exp(-((x-b)^2/(2*sigma^2))) regarding the definition of FWHM = sigma*2*sqrt(2*ln(2))
% therefore sigma = c1/sqrt(2) and –> FWHM = c1*2*sqrt(ln(2))
FWHM = (GaussFit.c2)*2*sqrt(log(2)); %[GHz]
fig2 = figure(‘Position’, [50 50 , 800, 600], ‘Color’, ‘white’);
plot(detuning,N_obs_2,’b’,’LineWidth’,2.5)
hold on
p2 = plot(GaussFit,’r:’);
p2.LineWidth = 2.7;
title([‘FWHM = ‘,num2str(FWHM),’GHz’])
xlabel(‘laser frequency [GHz]’)
ylabel(‘N_{obs}’)
set(gca,’FontSize’,18)
toc
%% SOLVING THE RATE EQNS %%
function [f_1_Sol,tspanOut,f_Sol] = solveLIFRateEqns(A_col,B_col,C_col,D_col,initialConds,tspan,time,p)
opts = odeset(‘RelTol’,1e-2,’AbsTol’,1e-4);
[tspanOut,f_Sol] = ode15s(@(t,f) makeDiff(t,f,time,A_col,B_col,C_col,D_col,p),tspan,initialConds,opts);
% THIS WILL SOLVE THE EQUAITON FOR EACH INITIAL CONDITION over time! so if there’s a
% 12×1 initial condition array, and a 100×1 time array, there will be
% 100×12 matrix output where each COLUMN is the solution (over time, i.e., the rows of said column) for each initial condition.
midpt = length(initialConds)/2;
f_0_Sol = f_Sol(:,1:midpt); % (n x p)
f_1_Sol = f_Sol(:,midpt + 1 :end); % (n x p)
function dfdt = makeDiff(t,f,time_interp,A_col,B_col,C_col,D_col,p)
% input:
% t = tspan array; independent variable (1 x 1)
% f = initial conditions array; dependent variable
%output:
% dfdt = an mx1 array used as the input for ode45
zero = (1:p);
one = ((p + 1):length(f));
% Interpolate the data sets (time,coeffMatrixX) at time t
aa = interp1(time_interp,A_col,t,’v5cubic’);
bb = interp1(time_interp,B_col,t,’v5cubic’);
cc = interp1(time_interp,C_col,t,’v5cubic’);
dd = interp1(time_interp,D_col,t,’v5cubic’);
dfdt = zeros(length(f),1);
dfdt(zero) = bb.*f(one) – aa.*f(zero);
dfdt(one) = cc.*f(zero) – dd.*f(one);
end
end
%% FITTING FUNCTION %%
function [fitresult, gof] = Gauss2(vLimits, N_obs_numerical)
%% Fit: ‘untitled fit 1’.
[xData, yData] = prepareCurveData( vLimits, N_obs_numerical );
% Set up fittype and options.
ft = fittype( ‘gauss2’ );
opts = fitoptions( ‘Method’, ‘NonlinearLeastSquares’ );
opts.Display = ‘Off’;
opts.Lower = [-Inf -Inf 0 -Inf -Inf 0];
opts.StartPoint = [418155525865.511 0.0071448985625 42.2244133622166 47891.96875 0 25];
% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
% Plot fit with data.
figure( ‘Name’, ‘untitled fit 1’ );
h = plot( fitresult, xData, yData );
legend( h, ‘N_obs_numerical vs. vLimits’, ‘untitled fit 1’, ‘Location’, ‘NorthEast’, ‘Interpreter’, ‘none’ );
% Label axes
xlabel( ‘vLimits’, ‘Interpreter’, ‘none’ );
ylabel( ‘N_obs_numerical’, ‘Interpreter’, ‘none’ );
grid on
endHello, I am not getting the correct answer when computing the ode for this set of rate equations. I believe the issue comes with the interpolation step because when I replace the coefficients with the non-interpolated values (A_col, B_col, C_col, D_col in place of aa, bb, cc, dd), I get the correct answer (a lorentzian). I am only using an array of ones as my "time dependent" part of the equation as I am using this as a test to see if I can get the correct answer before I make things more complicated.
I have tried all of the different interpolation schemes available (linear’, ‘nearest’, etc.). I have tried lowering the tolerance.
Any ideas?!
Thanks in advance!!
%clear;
%% ——————- CONSTANTS ——————- %%
tic
m_ArII = 6.6335e-26; %[kg]
e = 1.602e-19; %[C]
c = 299792458; %[m/s]
k_B = 1.380649e-23; %[J/K]
lambda_0 = 6.116616e-7;%[m]
freq_0 = c/lambda_0; %[Hz]
I_0 = 0.5e+9; %[W/m^2]
lineWidth = 1e+9; % [Hz]
A_1_0 = 0.212e+8; %[Hz]
A_1_2 = 0.759e+8; %[Hz]
B_1_0 = 121.9e+11; %[m^2/(J s)]
B_0_1 = 97.55e+11; %[m^2/(J s)]
solidAngle = sin(15*pi/180);
pulseWidth = 6e-9; % [s]
tau_0 = 1e-3;
tau_1 = 8.51e-9;
n_0 = 5e+13; %[1/m^3]
T_i = 0.1; %[eV]
p = 20; %number of intergrals
% —————————————————— %
%———- Limits for integration and diffeqs———-%
freqRange = 30e+9;
freqLimits = linspace(freq_0 – (freqRange*3),freq_0 + (freqRange*3),p); %[Hz]
centralLaserFreq = linspace((freq_0 – freqRange), (freq_0 + freqRange) ,p); %[Hz]
vLimits = linspace(-1e+4,1e+4,p)’; %[m/s]
tspan = [0 pulseWidth];
time_interp = linspace(0,pulseWidth,p)’;
%———————————————%
detuning = (centralLaserFreq – freq_0)./(1e+9);
% preallocated arrays
N_obs_2 = zeros(p,1);
f1intStore = zeros(p,p);
f_1_compare = zeros(p,p*p);
protoStore = zeros(p,p);
%——set up mesgrid for integration——-%
%rows %cols
[F , V] = meshgrid(freqLimits,vLimits);
%——————————————-%
%ion absorption spectrum
L = (A_1_0/pi)./((F – freq_0*(1 + (V/c)) ).^2 + (A_1_0)^2); %freqLimits is a (p x p) array
%———— INITIAL CONDITIONS ————%
f_0_initial = n_0*sqrt(m_ArII/(2*pi*e*T_i)).*exp(-(m_ArII*(vLimits).^2)./(2*e*T_i)); % (p x 1)
beta = (1e-7);
n_1 = beta*n_0;
f_1_initial = n_1*sqrt(m_ArII/(2*pi*e*T_i)).*exp(-(m_ArII*(vLimits).^2)./(2*e*T_i)); % (p x 1);
initialConds = [f_0_initial; f_1_initial]; %(2*p x 1)
%——————————————–%
%time component
theta = ones(p,1);% this is just an array of 1s for now, but will be a temporally dependent term in the future
%loop through different values of centralLaserFreq
for k = 1:p
%laser intensity spectrum
I_of_nu = I_0.*((lineWidth/(2*pi)))./((F – (centralLaserFreq(k) )).^2 + (lineWidth/2)^2);%(p x p)
%effective laser intensity
protoPhi = (1/(4*pi)).*I_of_nu.*L; %(p x p)
%integrate of frequency
protoPhiInt = trapz(freqLimits,protoPhi,2); %(p x 1)
protoStore(:,k) = protoPhiInt;
% (1 x p) (p x 1)
Phi_of_v_t = theta’.*protoPhiInt; %(p x p) % <— Now that the integration over nu is complete, we can bring in the theta term. Didn’t add it previously because it would’ve made things more complicated.
%example of Phi_of_v_t:
% [t1*phi1, t2*phi1, t3*phi1]
% [t1*phi2, t2*phi2, t3*phi2]
% [t1*phi3, t2*phi3, t3*phi3]
%Coeffitients to the rate equations
A = ((1/tau_0) + B_0_1.*Phi_of_v_t); %(p x p)
B = A_1_0 + B_1_0.*Phi_of_v_t; %(p x p)
C = B_0_1.*Phi_of_v_t; %(p x p)
D = ((1/tau_1) + B_1_0.*Phi_of_v_t); %(p x p)
%loop through columns of the coefficients to make things faster than going
%through each element individually
f1intStore = zeros(p,p);
for col = 1:p
[f_1_Sol,tspanOut,f_Sol] = solveLIFRateEqns(A(:,col),B(:,col),C(:,col),D(:,col),initialConds,tspan,time_interp,p);
f_1_timeInt = trapz(tspanOut,f_1_Sol,1); %integrate over time (1 x p)
f1intStore(:,col) = f_1_timeInt; % each column of f1solStore is the f_1 for a different time
end
N_obs_int = (solidAngle/(4*pi)).*A_1_2.*f1intStore(:,1); %(p x 1)
fxintAll2 = trapz(vLimits,N_obs_int); %intgrate over veocity
N_obs_2(k,1) = fxintAll2; %store value
disp(k)
end
%% PLOTTING
[GaussFit, GaussGof] = Gauss2(detuning’, N_obs_2);
% MATLAB’s "gauss1" fitting defines the Gaussian as a1*exp(-((x-b1)/c1)^2) and gives c1 as a result of fitting.
% I’m using a Gaussian defined as a*exp(-((x-b)^2/(2*sigma^2))) regarding the definition of FWHM = sigma*2*sqrt(2*ln(2))
% therefore sigma = c1/sqrt(2) and –> FWHM = c1*2*sqrt(ln(2))
FWHM = (GaussFit.c2)*2*sqrt(log(2)); %[GHz]
fig2 = figure(‘Position’, [50 50 , 800, 600], ‘Color’, ‘white’);
plot(detuning,N_obs_2,’b’,’LineWidth’,2.5)
hold on
p2 = plot(GaussFit,’r:’);
p2.LineWidth = 2.7;
title([‘FWHM = ‘,num2str(FWHM),’GHz’])
xlabel(‘laser frequency [GHz]’)
ylabel(‘N_{obs}’)
set(gca,’FontSize’,18)
toc
%% SOLVING THE RATE EQNS %%
function [f_1_Sol,tspanOut,f_Sol] = solveLIFRateEqns(A_col,B_col,C_col,D_col,initialConds,tspan,time,p)
opts = odeset(‘RelTol’,1e-2,’AbsTol’,1e-4);
[tspanOut,f_Sol] = ode15s(@(t,f) makeDiff(t,f,time,A_col,B_col,C_col,D_col,p),tspan,initialConds,opts);
% THIS WILL SOLVE THE EQUAITON FOR EACH INITIAL CONDITION over time! so if there’s a
% 12×1 initial condition array, and a 100×1 time array, there will be
% 100×12 matrix output where each COLUMN is the solution (over time, i.e., the rows of said column) for each initial condition.
midpt = length(initialConds)/2;
f_0_Sol = f_Sol(:,1:midpt); % (n x p)
f_1_Sol = f_Sol(:,midpt + 1 :end); % (n x p)
function dfdt = makeDiff(t,f,time_interp,A_col,B_col,C_col,D_col,p)
% input:
% t = tspan array; independent variable (1 x 1)
% f = initial conditions array; dependent variable
%output:
% dfdt = an mx1 array used as the input for ode45
zero = (1:p);
one = ((p + 1):length(f));
% Interpolate the data sets (time,coeffMatrixX) at time t
aa = interp1(time_interp,A_col,t,’v5cubic’);
bb = interp1(time_interp,B_col,t,’v5cubic’);
cc = interp1(time_interp,C_col,t,’v5cubic’);
dd = interp1(time_interp,D_col,t,’v5cubic’);
dfdt = zeros(length(f),1);
dfdt(zero) = bb.*f(one) – aa.*f(zero);
dfdt(one) = cc.*f(zero) – dd.*f(one);
end
end
%% FITTING FUNCTION %%
function [fitresult, gof] = Gauss2(vLimits, N_obs_numerical)
%% Fit: ‘untitled fit 1’.
[xData, yData] = prepareCurveData( vLimits, N_obs_numerical );
% Set up fittype and options.
ft = fittype( ‘gauss2’ );
opts = fitoptions( ‘Method’, ‘NonlinearLeastSquares’ );
opts.Display = ‘Off’;
opts.Lower = [-Inf -Inf 0 -Inf -Inf 0];
opts.StartPoint = [418155525865.511 0.0071448985625 42.2244133622166 47891.96875 0 25];
% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
% Plot fit with data.
figure( ‘Name’, ‘untitled fit 1’ );
h = plot( fitresult, xData, yData );
legend( h, ‘N_obs_numerical vs. vLimits’, ‘untitled fit 1’, ‘Location’, ‘NorthEast’, ‘Interpreter’, ‘none’ );
% Label axes
xlabel( ‘vLimits’, ‘Interpreter’, ‘none’ );
ylabel( ‘N_obs_numerical’, ‘Interpreter’, ‘none’ );
grid on
end Hello, I am not getting the correct answer when computing the ode for this set of rate equations. I believe the issue comes with the interpolation step because when I replace the coefficients with the non-interpolated values (A_col, B_col, C_col, D_col in place of aa, bb, cc, dd), I get the correct answer (a lorentzian). I am only using an array of ones as my "time dependent" part of the equation as I am using this as a test to see if I can get the correct answer before I make things more complicated.
I have tried all of the different interpolation schemes available (linear’, ‘nearest’, etc.). I have tried lowering the tolerance.
Any ideas?!
Thanks in advance!!
%clear;
%% ——————- CONSTANTS ——————- %%
tic
m_ArII = 6.6335e-26; %[kg]
e = 1.602e-19; %[C]
c = 299792458; %[m/s]
k_B = 1.380649e-23; %[J/K]
lambda_0 = 6.116616e-7;%[m]
freq_0 = c/lambda_0; %[Hz]
I_0 = 0.5e+9; %[W/m^2]
lineWidth = 1e+9; % [Hz]
A_1_0 = 0.212e+8; %[Hz]
A_1_2 = 0.759e+8; %[Hz]
B_1_0 = 121.9e+11; %[m^2/(J s)]
B_0_1 = 97.55e+11; %[m^2/(J s)]
solidAngle = sin(15*pi/180);
pulseWidth = 6e-9; % [s]
tau_0 = 1e-3;
tau_1 = 8.51e-9;
n_0 = 5e+13; %[1/m^3]
T_i = 0.1; %[eV]
p = 20; %number of intergrals
% —————————————————— %
%———- Limits for integration and diffeqs———-%
freqRange = 30e+9;
freqLimits = linspace(freq_0 – (freqRange*3),freq_0 + (freqRange*3),p); %[Hz]
centralLaserFreq = linspace((freq_0 – freqRange), (freq_0 + freqRange) ,p); %[Hz]
vLimits = linspace(-1e+4,1e+4,p)’; %[m/s]
tspan = [0 pulseWidth];
time_interp = linspace(0,pulseWidth,p)’;
%———————————————%
detuning = (centralLaserFreq – freq_0)./(1e+9);
% preallocated arrays
N_obs_2 = zeros(p,1);
f1intStore = zeros(p,p);
f_1_compare = zeros(p,p*p);
protoStore = zeros(p,p);
%——set up mesgrid for integration——-%
%rows %cols
[F , V] = meshgrid(freqLimits,vLimits);
%——————————————-%
%ion absorption spectrum
L = (A_1_0/pi)./((F – freq_0*(1 + (V/c)) ).^2 + (A_1_0)^2); %freqLimits is a (p x p) array
%———— INITIAL CONDITIONS ————%
f_0_initial = n_0*sqrt(m_ArII/(2*pi*e*T_i)).*exp(-(m_ArII*(vLimits).^2)./(2*e*T_i)); % (p x 1)
beta = (1e-7);
n_1 = beta*n_0;
f_1_initial = n_1*sqrt(m_ArII/(2*pi*e*T_i)).*exp(-(m_ArII*(vLimits).^2)./(2*e*T_i)); % (p x 1);
initialConds = [f_0_initial; f_1_initial]; %(2*p x 1)
%——————————————–%
%time component
theta = ones(p,1);% this is just an array of 1s for now, but will be a temporally dependent term in the future
%loop through different values of centralLaserFreq
for k = 1:p
%laser intensity spectrum
I_of_nu = I_0.*((lineWidth/(2*pi)))./((F – (centralLaserFreq(k) )).^2 + (lineWidth/2)^2);%(p x p)
%effective laser intensity
protoPhi = (1/(4*pi)).*I_of_nu.*L; %(p x p)
%integrate of frequency
protoPhiInt = trapz(freqLimits,protoPhi,2); %(p x 1)
protoStore(:,k) = protoPhiInt;
% (1 x p) (p x 1)
Phi_of_v_t = theta’.*protoPhiInt; %(p x p) % <— Now that the integration over nu is complete, we can bring in the theta term. Didn’t add it previously because it would’ve made things more complicated.
%example of Phi_of_v_t:
% [t1*phi1, t2*phi1, t3*phi1]
% [t1*phi2, t2*phi2, t3*phi2]
% [t1*phi3, t2*phi3, t3*phi3]
%Coeffitients to the rate equations
A = ((1/tau_0) + B_0_1.*Phi_of_v_t); %(p x p)
B = A_1_0 + B_1_0.*Phi_of_v_t; %(p x p)
C = B_0_1.*Phi_of_v_t; %(p x p)
D = ((1/tau_1) + B_1_0.*Phi_of_v_t); %(p x p)
%loop through columns of the coefficients to make things faster than going
%through each element individually
f1intStore = zeros(p,p);
for col = 1:p
[f_1_Sol,tspanOut,f_Sol] = solveLIFRateEqns(A(:,col),B(:,col),C(:,col),D(:,col),initialConds,tspan,time_interp,p);
f_1_timeInt = trapz(tspanOut,f_1_Sol,1); %integrate over time (1 x p)
f1intStore(:,col) = f_1_timeInt; % each column of f1solStore is the f_1 for a different time
end
N_obs_int = (solidAngle/(4*pi)).*A_1_2.*f1intStore(:,1); %(p x 1)
fxintAll2 = trapz(vLimits,N_obs_int); %intgrate over veocity
N_obs_2(k,1) = fxintAll2; %store value
disp(k)
end
%% PLOTTING
[GaussFit, GaussGof] = Gauss2(detuning’, N_obs_2);
% MATLAB’s "gauss1" fitting defines the Gaussian as a1*exp(-((x-b1)/c1)^2) and gives c1 as a result of fitting.
% I’m using a Gaussian defined as a*exp(-((x-b)^2/(2*sigma^2))) regarding the definition of FWHM = sigma*2*sqrt(2*ln(2))
% therefore sigma = c1/sqrt(2) and –> FWHM = c1*2*sqrt(ln(2))
FWHM = (GaussFit.c2)*2*sqrt(log(2)); %[GHz]
fig2 = figure(‘Position’, [50 50 , 800, 600], ‘Color’, ‘white’);
plot(detuning,N_obs_2,’b’,’LineWidth’,2.5)
hold on
p2 = plot(GaussFit,’r:’);
p2.LineWidth = 2.7;
title([‘FWHM = ‘,num2str(FWHM),’GHz’])
xlabel(‘laser frequency [GHz]’)
ylabel(‘N_{obs}’)
set(gca,’FontSize’,18)
toc
%% SOLVING THE RATE EQNS %%
function [f_1_Sol,tspanOut,f_Sol] = solveLIFRateEqns(A_col,B_col,C_col,D_col,initialConds,tspan,time,p)
opts = odeset(‘RelTol’,1e-2,’AbsTol’,1e-4);
[tspanOut,f_Sol] = ode15s(@(t,f) makeDiff(t,f,time,A_col,B_col,C_col,D_col,p),tspan,initialConds,opts);
% THIS WILL SOLVE THE EQUAITON FOR EACH INITIAL CONDITION over time! so if there’s a
% 12×1 initial condition array, and a 100×1 time array, there will be
% 100×12 matrix output where each COLUMN is the solution (over time, i.e., the rows of said column) for each initial condition.
midpt = length(initialConds)/2;
f_0_Sol = f_Sol(:,1:midpt); % (n x p)
f_1_Sol = f_Sol(:,midpt + 1 :end); % (n x p)
function dfdt = makeDiff(t,f,time_interp,A_col,B_col,C_col,D_col,p)
% input:
% t = tspan array; independent variable (1 x 1)
% f = initial conditions array; dependent variable
%output:
% dfdt = an mx1 array used as the input for ode45
zero = (1:p);
one = ((p + 1):length(f));
% Interpolate the data sets (time,coeffMatrixX) at time t
aa = interp1(time_interp,A_col,t,’v5cubic’);
bb = interp1(time_interp,B_col,t,’v5cubic’);
cc = interp1(time_interp,C_col,t,’v5cubic’);
dd = interp1(time_interp,D_col,t,’v5cubic’);
dfdt = zeros(length(f),1);
dfdt(zero) = bb.*f(one) – aa.*f(zero);
dfdt(one) = cc.*f(zero) – dd.*f(one);
end
end
%% FITTING FUNCTION %%
function [fitresult, gof] = Gauss2(vLimits, N_obs_numerical)
%% Fit: ‘untitled fit 1’.
[xData, yData] = prepareCurveData( vLimits, N_obs_numerical );
% Set up fittype and options.
ft = fittype( ‘gauss2’ );
opts = fitoptions( ‘Method’, ‘NonlinearLeastSquares’ );
opts.Display = ‘Off’;
opts.Lower = [-Inf -Inf 0 -Inf -Inf 0];
opts.StartPoint = [418155525865.511 0.0071448985625 42.2244133622166 47891.96875 0 25];
% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );
% Plot fit with data.
figure( ‘Name’, ‘untitled fit 1’ );
h = plot( fitresult, xData, yData );
legend( h, ‘N_obs_numerical vs. vLimits’, ‘untitled fit 1’, ‘Location’, ‘NorthEast’, ‘Interpreter’, ‘none’ );
% Label axes
xlabel( ‘vLimits’, ‘Interpreter’, ‘none’ );
ylabel( ‘N_obs_numerical’, ‘Interpreter’, ‘none’ );
grid on
end ode45s, differential equations, time-varying coefficients, rate equations, ode15s MATLAB Answers — New Questions
2021b & RHEL8
Does RHEL8 support 2021b? We are having an issue getting 2021b to install on RHEL8 with the following errors:
We are mounting the .iso using mount -o loop <iso file name> /mnt/iso and we are getting this error:
/mnt/iso/bin/glnxa64/MathWorksProductInstaller: error while loading shared libraries: libmwfoundation_log.so: cannot open shared object file: No such file or directory
If you do a search for the libmwfoundation_log.so you can find it on the mounted CD so not sure why it thinks it’s not there.Does RHEL8 support 2021b? We are having an issue getting 2021b to install on RHEL8 with the following errors:
We are mounting the .iso using mount -o loop <iso file name> /mnt/iso and we are getting this error:
/mnt/iso/bin/glnxa64/MathWorksProductInstaller: error while loading shared libraries: libmwfoundation_log.so: cannot open shared object file: No such file or directory
If you do a search for the libmwfoundation_log.so you can find it on the mounted CD so not sure why it thinks it’s not there. Does RHEL8 support 2021b? We are having an issue getting 2021b to install on RHEL8 with the following errors:
We are mounting the .iso using mount -o loop <iso file name> /mnt/iso and we are getting this error:
/mnt/iso/bin/glnxa64/MathWorksProductInstaller: error while loading shared libraries: libmwfoundation_log.so: cannot open shared object file: No such file or directory
If you do a search for the libmwfoundation_log.so you can find it on the mounted CD so not sure why it thinks it’s not there. rhel8, 2021b, matlab, installation MATLAB Answers — New Questions
Non-homogenous spatial and temporal grids
Good day family. I am solving a system of 2 PDEs using the pdepe package (The algorithm is in the matlab documentation: https://www.mathworks.com/help/matlab/math/solve-system-of-pdes.html ). The parameters of the system depend on some observations (data) and the spatial axis as well as shown in my code below:
% beggining of the code
function non_constant_spatial_and_temporal_axis
% pdebc function — Neumann BCs
function [pl, ql, pr, qr] = pdebc(~, ~, ~, ~, ~)
pl = [0; 0];
ql = [1; 1];
pr = [0; 0];
qr = [1;1];
end
% pdeic function
function u0 = pdeic(~)
u0 = [10;5];
end
% pdefun function
function [c, f, s] = pdefun(x, ~, u, dudx)
D1 = 0.024;
D2 = 0.0170;
c = [1; 1];
f = [D1; D2].*dudx;
t1=[1,0,3,4,5,3,2,3,2,3,3,2,3,4,5,6,7,5,4,3,3,2,2,1,3,2,4,5]; % This data is the observations which represent the temporal grid
x_interpolate = linspace(0, 10, length(t1)); % Interpolation of the data so that it can match with the spatial grid
t1_interpolated = interp1(x_interpolate,t1,x,’linear’);
alpha1 = 0.5*(1-sin(2*t1_interpolated)); % The parameters alpha1 and alpha2 depends on the observations in t1
% I want to define alpha1 in such a way that it is a function of bth
% time and space x such that alpha1(t,x)=
% 0.5*(1-sin(2*t1_interpolated))*(1-sin(2*x)) so that I can visualize
% the oscillations on both the time and the space axis at the same
% time.
alpha2 = 0.5*(1-sin(2*t1_interpolated)); % The parameter alpha2 depends on the observations in t1
% alpha2(t,x)=0.5*(1-sin(2*t1_interpolated))*(1-sin(2*x));
s = [alpha1.*u(2).*(u(2).^2 -1); -alpha2.* u(1).*(u(2)-1)];
end
% Main script for running the code
clear all; close all; clc;
m = 0;
t = linspace(0, 5, 50); % time variable
x = linspace(0, 10,50); % space varable
sol = pdepe(m, @pdefun, @pdeic, @pdebc, x, t);
% Extract the solutions
u1 = sol(:, :, 1);
u2 = sol(:, :, 2);
% Plots
figure(1)
surf(x, t, u1)
title(‘u_1(x,t)’)
xlabel(‘Time t’)
ylabel(‘Postion x’)
figure(2)
surf(x, t, u2)
title(‘u_2(x,t)’)
xlabel(‘Time t’)
ylabel(‘Position x’)
end
%end of code
The code is running well but its not giving me oscillations on the spatial axis, like what I am getting on the temporal axis. I wish to be able to observe the spatial variation using a sinusoidal function like the sin function and I also view the oscillations on the spatial and temporal axis. Note that the oscillations on the temporal axis are due to the data while the ones on the spatial axis are due to the use of the sinusoidal function. May you kindly assist me resolve this. I have tried substituting alpha1 with alpha1=0.5*(1-sin(2*t1_interpolated))*(1-sin(2*x)); but then I observe that the term (1-sin(2*x)) reduces the value of alpha1 and the effect is only felt on the temporal axis and not in the spatial axis. What can I do to have the effect felt on the spatial axis as well?Good day family. I am solving a system of 2 PDEs using the pdepe package (The algorithm is in the matlab documentation: https://www.mathworks.com/help/matlab/math/solve-system-of-pdes.html ). The parameters of the system depend on some observations (data) and the spatial axis as well as shown in my code below:
% beggining of the code
function non_constant_spatial_and_temporal_axis
% pdebc function — Neumann BCs
function [pl, ql, pr, qr] = pdebc(~, ~, ~, ~, ~)
pl = [0; 0];
ql = [1; 1];
pr = [0; 0];
qr = [1;1];
end
% pdeic function
function u0 = pdeic(~)
u0 = [10;5];
end
% pdefun function
function [c, f, s] = pdefun(x, ~, u, dudx)
D1 = 0.024;
D2 = 0.0170;
c = [1; 1];
f = [D1; D2].*dudx;
t1=[1,0,3,4,5,3,2,3,2,3,3,2,3,4,5,6,7,5,4,3,3,2,2,1,3,2,4,5]; % This data is the observations which represent the temporal grid
x_interpolate = linspace(0, 10, length(t1)); % Interpolation of the data so that it can match with the spatial grid
t1_interpolated = interp1(x_interpolate,t1,x,’linear’);
alpha1 = 0.5*(1-sin(2*t1_interpolated)); % The parameters alpha1 and alpha2 depends on the observations in t1
% I want to define alpha1 in such a way that it is a function of bth
% time and space x such that alpha1(t,x)=
% 0.5*(1-sin(2*t1_interpolated))*(1-sin(2*x)) so that I can visualize
% the oscillations on both the time and the space axis at the same
% time.
alpha2 = 0.5*(1-sin(2*t1_interpolated)); % The parameter alpha2 depends on the observations in t1
% alpha2(t,x)=0.5*(1-sin(2*t1_interpolated))*(1-sin(2*x));
s = [alpha1.*u(2).*(u(2).^2 -1); -alpha2.* u(1).*(u(2)-1)];
end
% Main script for running the code
clear all; close all; clc;
m = 0;
t = linspace(0, 5, 50); % time variable
x = linspace(0, 10,50); % space varable
sol = pdepe(m, @pdefun, @pdeic, @pdebc, x, t);
% Extract the solutions
u1 = sol(:, :, 1);
u2 = sol(:, :, 2);
% Plots
figure(1)
surf(x, t, u1)
title(‘u_1(x,t)’)
xlabel(‘Time t’)
ylabel(‘Postion x’)
figure(2)
surf(x, t, u2)
title(‘u_2(x,t)’)
xlabel(‘Time t’)
ylabel(‘Position x’)
end
%end of code
The code is running well but its not giving me oscillations on the spatial axis, like what I am getting on the temporal axis. I wish to be able to observe the spatial variation using a sinusoidal function like the sin function and I also view the oscillations on the spatial and temporal axis. Note that the oscillations on the temporal axis are due to the data while the ones on the spatial axis are due to the use of the sinusoidal function. May you kindly assist me resolve this. I have tried substituting alpha1 with alpha1=0.5*(1-sin(2*t1_interpolated))*(1-sin(2*x)); but then I observe that the term (1-sin(2*x)) reduces the value of alpha1 and the effect is only felt on the temporal axis and not in the spatial axis. What can I do to have the effect felt on the spatial axis as well? Good day family. I am solving a system of 2 PDEs using the pdepe package (The algorithm is in the matlab documentation: https://www.mathworks.com/help/matlab/math/solve-system-of-pdes.html ). The parameters of the system depend on some observations (data) and the spatial axis as well as shown in my code below:
% beggining of the code
function non_constant_spatial_and_temporal_axis
% pdebc function — Neumann BCs
function [pl, ql, pr, qr] = pdebc(~, ~, ~, ~, ~)
pl = [0; 0];
ql = [1; 1];
pr = [0; 0];
qr = [1;1];
end
% pdeic function
function u0 = pdeic(~)
u0 = [10;5];
end
% pdefun function
function [c, f, s] = pdefun(x, ~, u, dudx)
D1 = 0.024;
D2 = 0.0170;
c = [1; 1];
f = [D1; D2].*dudx;
t1=[1,0,3,4,5,3,2,3,2,3,3,2,3,4,5,6,7,5,4,3,3,2,2,1,3,2,4,5]; % This data is the observations which represent the temporal grid
x_interpolate = linspace(0, 10, length(t1)); % Interpolation of the data so that it can match with the spatial grid
t1_interpolated = interp1(x_interpolate,t1,x,’linear’);
alpha1 = 0.5*(1-sin(2*t1_interpolated)); % The parameters alpha1 and alpha2 depends on the observations in t1
% I want to define alpha1 in such a way that it is a function of bth
% time and space x such that alpha1(t,x)=
% 0.5*(1-sin(2*t1_interpolated))*(1-sin(2*x)) so that I can visualize
% the oscillations on both the time and the space axis at the same
% time.
alpha2 = 0.5*(1-sin(2*t1_interpolated)); % The parameter alpha2 depends on the observations in t1
% alpha2(t,x)=0.5*(1-sin(2*t1_interpolated))*(1-sin(2*x));
s = [alpha1.*u(2).*(u(2).^2 -1); -alpha2.* u(1).*(u(2)-1)];
end
% Main script for running the code
clear all; close all; clc;
m = 0;
t = linspace(0, 5, 50); % time variable
x = linspace(0, 10,50); % space varable
sol = pdepe(m, @pdefun, @pdeic, @pdebc, x, t);
% Extract the solutions
u1 = sol(:, :, 1);
u2 = sol(:, :, 2);
% Plots
figure(1)
surf(x, t, u1)
title(‘u_1(x,t)’)
xlabel(‘Time t’)
ylabel(‘Postion x’)
figure(2)
surf(x, t, u2)
title(‘u_2(x,t)’)
xlabel(‘Time t’)
ylabel(‘Position x’)
end
%end of code
The code is running well but its not giving me oscillations on the spatial axis, like what I am getting on the temporal axis. I wish to be able to observe the spatial variation using a sinusoidal function like the sin function and I also view the oscillations on the spatial and temporal axis. Note that the oscillations on the temporal axis are due to the data while the ones on the spatial axis are due to the use of the sinusoidal function. May you kindly assist me resolve this. I have tried substituting alpha1 with alpha1=0.5*(1-sin(2*t1_interpolated))*(1-sin(2*x)); but then I observe that the term (1-sin(2*x)) reduces the value of alpha1 and the effect is only felt on the temporal axis and not in the spatial axis. What can I do to have the effect felt on the spatial axis as well? non-homogeneous spatial axis MATLAB Answers — New Questions
Can I simulate more than one stl geometry/ an assmebly with the ‘createpde’ function?
Hello, so I want to run a thermal (transient) thermal simulation on a CAD model.
I have imported the CAD model, which consists of distinct bodies that are in both thermal contact and physical contact (like a cuboid, inside another cubiod). I want to use the createpde fucntion to carryout the simulation.
This is an example of my code:
xxChamber = createpde(‘thermal’, ‘transient’);
importGeometry(xxChamber, ‘Outermost component v.stl’)
pdegplot(xxChamber, "FaceLabels","on")
Now this works quite alright, but once I begin to add the second chamber, by typing in
importGeometry (hydrolysisChamber, "middle v2.stl")
The code gives an error saying ‘ Assemblies of more than one geometric model not supported.
What can I do to solve this problem?Hello, so I want to run a thermal (transient) thermal simulation on a CAD model.
I have imported the CAD model, which consists of distinct bodies that are in both thermal contact and physical contact (like a cuboid, inside another cubiod). I want to use the createpde fucntion to carryout the simulation.
This is an example of my code:
xxChamber = createpde(‘thermal’, ‘transient’);
importGeometry(xxChamber, ‘Outermost component v.stl’)
pdegplot(xxChamber, "FaceLabels","on")
Now this works quite alright, but once I begin to add the second chamber, by typing in
importGeometry (hydrolysisChamber, "middle v2.stl")
The code gives an error saying ‘ Assemblies of more than one geometric model not supported.
What can I do to solve this problem? Hello, so I want to run a thermal (transient) thermal simulation on a CAD model.
I have imported the CAD model, which consists of distinct bodies that are in both thermal contact and physical contact (like a cuboid, inside another cubiod). I want to use the createpde fucntion to carryout the simulation.
This is an example of my code:
xxChamber = createpde(‘thermal’, ‘transient’);
importGeometry(xxChamber, ‘Outermost component v.stl’)
pdegplot(xxChamber, "FaceLabels","on")
Now this works quite alright, but once I begin to add the second chamber, by typing in
importGeometry (hydrolysisChamber, "middle v2.stl")
The code gives an error saying ‘ Assemblies of more than one geometric model not supported.
What can I do to solve this problem? pde, simulation, model MATLAB Answers — New Questions
Transfer function derviation for a system
I plan to heat 10L of Palm Oil to 60C from room temperature and maintain that temperature for the transesterification of biodiesel. To maintain the temperature I am developing a PID controller with a 1000W immersion heater but I need to find a way to obtain the PID constants. So far I have developed this equation to model:
Where P is input power,
How to
derive the system transfer function from matlab
can I generate similar modelled equations as above for any adjustments to the system for example if I add a stirrer?I plan to heat 10L of Palm Oil to 60C from room temperature and maintain that temperature for the transesterification of biodiesel. To maintain the temperature I am developing a PID controller with a 1000W immersion heater but I need to find a way to obtain the PID constants. So far I have developed this equation to model:
Where P is input power,
How to
derive the system transfer function from matlab
can I generate similar modelled equations as above for any adjustments to the system for example if I add a stirrer? I plan to heat 10L of Palm Oil to 60C from room temperature and maintain that temperature for the transesterification of biodiesel. To maintain the temperature I am developing a PID controller with a 1000W immersion heater but I need to find a way to obtain the PID constants. So far I have developed this equation to model:
Where P is input power,
How to
derive the system transfer function from matlab
can I generate similar modelled equations as above for any adjustments to the system for example if I add a stirrer? transfer, thermo MATLAB Answers — New Questions
add / interpolate arrays of different length
Hello community,
i have two array’s tsol and ysol. Each array has the length of 9;
I have also an error estimate array errorEst. The error estimate has only 4 elements.
The four elements have a "time" length of 0.25.
The errorEst
tsol = linspace(0,1,9);
ysol = sin(tsol*pi);
errorEst = rand(1,4); %each bar has a length of 0.25;
bar(errorEst);
plot([0 0.25 0.5 0.75; 0.25 0.5 0.75 1],[errorEst; errorEst]/0.25,’k’)
title(‘gradient of errorEst ‘)
xlabel(‘t’)
ylabel(‘d errorEst / dt’)
% who could i integrate this gradient of errorEst w.r.t the time tSol and
% add to ysol?
% the first part should be a continous (function/plot)
% is there an elegant solution, which is parallelizable?
the solution should be plottet:
%yp = ysol + errorEstLin;
%ym= ysol -errorEstLin;Hello community,
i have two array’s tsol and ysol. Each array has the length of 9;
I have also an error estimate array errorEst. The error estimate has only 4 elements.
The four elements have a "time" length of 0.25.
The errorEst
tsol = linspace(0,1,9);
ysol = sin(tsol*pi);
errorEst = rand(1,4); %each bar has a length of 0.25;
bar(errorEst);
plot([0 0.25 0.5 0.75; 0.25 0.5 0.75 1],[errorEst; errorEst]/0.25,’k’)
title(‘gradient of errorEst ‘)
xlabel(‘t’)
ylabel(‘d errorEst / dt’)
% who could i integrate this gradient of errorEst w.r.t the time tSol and
% add to ysol?
% the first part should be a continous (function/plot)
% is there an elegant solution, which is parallelizable?
the solution should be plottet:
%yp = ysol + errorEstLin;
%ym= ysol -errorEstLin; Hello community,
i have two array’s tsol and ysol. Each array has the length of 9;
I have also an error estimate array errorEst. The error estimate has only 4 elements.
The four elements have a "time" length of 0.25.
The errorEst
tsol = linspace(0,1,9);
ysol = sin(tsol*pi);
errorEst = rand(1,4); %each bar has a length of 0.25;
bar(errorEst);
plot([0 0.25 0.5 0.75; 0.25 0.5 0.75 1],[errorEst; errorEst]/0.25,’k’)
title(‘gradient of errorEst ‘)
xlabel(‘t’)
ylabel(‘d errorEst / dt’)
% who could i integrate this gradient of errorEst w.r.t the time tSol and
% add to ysol?
% the first part should be a continous (function/plot)
% is there an elegant solution, which is parallelizable?
the solution should be plottet:
%yp = ysol + errorEstLin;
%ym= ysol -errorEstLin; interpolate MATLAB Answers — New Questions
How do I resolve the error “java.lang.UnsatisfiedLinkError: nativemvm.dll already loaded in another classloader” when running a Java Engine program?
I have a Java source code that calls MATLAB via MATLAB Engine API. I compiled a .WAR file from the Java source code on Windows 7. Then I deployed the .WAR file to Windows Server 2016 using JBoss EAP, and got the following error:
ERROR: java.lang.UnsatisfiedLinkError: Native Library C:Program FilesMATLABR2019abinwin64nativemvm.dll already loaded in another classloader
I have MATLAB R2019a installed on both Windows 7 and Windows Sever 2016. Before deploying the .WAR file, I added "C:Program FilesMATLABR2019abinwin64" to the PATH variable on Windows Server 2016 and added “C:Program FilesMATLABR2019aexternenginesjavajarengine.jar” to the “classpath” in JBoss.
How do I resolve this "nativemvm.dll" error?I have a Java source code that calls MATLAB via MATLAB Engine API. I compiled a .WAR file from the Java source code on Windows 7. Then I deployed the .WAR file to Windows Server 2016 using JBoss EAP, and got the following error:
ERROR: java.lang.UnsatisfiedLinkError: Native Library C:Program FilesMATLABR2019abinwin64nativemvm.dll already loaded in another classloader
I have MATLAB R2019a installed on both Windows 7 and Windows Sever 2016. Before deploying the .WAR file, I added "C:Program FilesMATLABR2019abinwin64" to the PATH variable on Windows Server 2016 and added “C:Program FilesMATLABR2019aexternenginesjavajarengine.jar” to the “classpath” in JBoss.
How do I resolve this "nativemvm.dll" error? I have a Java source code that calls MATLAB via MATLAB Engine API. I compiled a .WAR file from the Java source code on Windows 7. Then I deployed the .WAR file to Windows Server 2016 using JBoss EAP, and got the following error:
ERROR: java.lang.UnsatisfiedLinkError: Native Library C:Program FilesMATLABR2019abinwin64nativemvm.dll already loaded in another classloader
I have MATLAB R2019a installed on both Windows 7 and Windows Sever 2016. Before deploying the .WAR file, I added "C:Program FilesMATLABR2019abinwin64" to the PATH variable on Windows Server 2016 and added “C:Program FilesMATLABR2019aexternenginesjavajarengine.jar” to the “classpath” in JBoss.
How do I resolve this "nativemvm.dll" error? MATLAB Answers — New Questions
How to choose relative error bound when using reducespec on sparse LTI model
Hello everyone,
I have a sparse LTI model, named "sys", that I want to reduce with reducespec and balanced truncation.
The final model should have a maximum relative error of e.g. 30%.
I have tried the following:
R = reducespec(sys, "balanced");
R.Options.Goal = "relative";
sysR = getrom(R, MaxError=0.3, Method=’truncate’);
However, in contrast to the full model, the MOR specification object of the sparse model does not have the "Goal" field under "R.Options", which appears to mean that it cannot be switched to relative error bound.
Unrecognized property ‘Goal’ for class ‘mor.SparseBalancedTruncationOptions’.
Is there another way to switch to the relative error bound, or is there a fundamental underlying reason why this is not applicable to sparse models?
Thank you in advance!Hello everyone,
I have a sparse LTI model, named "sys", that I want to reduce with reducespec and balanced truncation.
The final model should have a maximum relative error of e.g. 30%.
I have tried the following:
R = reducespec(sys, "balanced");
R.Options.Goal = "relative";
sysR = getrom(R, MaxError=0.3, Method=’truncate’);
However, in contrast to the full model, the MOR specification object of the sparse model does not have the "Goal" field under "R.Options", which appears to mean that it cannot be switched to relative error bound.
Unrecognized property ‘Goal’ for class ‘mor.SparseBalancedTruncationOptions’.
Is there another way to switch to the relative error bound, or is there a fundamental underlying reason why this is not applicable to sparse models?
Thank you in advance! Hello everyone,
I have a sparse LTI model, named "sys", that I want to reduce with reducespec and balanced truncation.
The final model should have a maximum relative error of e.g. 30%.
I have tried the following:
R = reducespec(sys, "balanced");
R.Options.Goal = "relative";
sysR = getrom(R, MaxError=0.3, Method=’truncate’);
However, in contrast to the full model, the MOR specification object of the sparse model does not have the "Goal" field under "R.Options", which appears to mean that it cannot be switched to relative error bound.
Unrecognized property ‘Goal’ for class ‘mor.SparseBalancedTruncationOptions’.
Is there another way to switch to the relative error bound, or is there a fundamental underlying reason why this is not applicable to sparse models?
Thank you in advance! reducespec, sparse lti, relative error bound MATLAB Answers — New Questions
Troubles With Image Resizing
Hello! I have a slice of a .nii file (essentially a 3D array) I’d like to print out. However, the dimensions of a voxel (a 3D pixel, so basically a cube) aren’t exactly 1mm x 1mm x 1mm. They’re around 1.2mm x 1.1mm x 1.1mm.
I took a slice of the 3D array, so it’s just a 2D image. So, the dimensions of each pixel are ~1.2mm x 1.1mm. I’d like to resize the image such that one pixel is equal to one element of the array is equal to one mm. In theory, it doesn’t seem too difficult. But I just can’t figure this out.
Here’s my code so far:
clear
clc
close all
% Load NIfTI image
info = niftiinfo(‘index_image.nii’);
brain = niftiread(info);
% Extract a slice from the NIfTI image
slice = brain(:, :, 128);
% Create figure
hFig = figure;
% Display the slice with pixel-to-mm scaling
imagesc(slice);
axis off;
axis equal;
colormap(gray)
% Set figure size so that 1 pixel corresponds to 1mm on the printed figure
set(hFig, ‘Units’, ‘centimeters’, ‘Position’, [0 0 24 17.6]);
movegui(hFig, ‘center’);
% Set properties to control the output size
set(hFig, ‘PaperUnits’, ‘centimeters’);
set(hFig, ‘PaperPosition’, [0 0 24 17.6]);
set(hFig, ‘PaperOrientation’, ‘landscape’);
% Export to PDF and open file
print(hFig, ‘-dpdf’, ‘-r0’, ‘out.pdf’);
open(‘out.pdf’);
The reason why use 24 and 17.6 is because my .nii file is 176 x 240 x 256 elements. I’d honestly like the code to work regardless of the size of the array, but even hard coding it into MATLAB isn’t working. The code is compiling just fine, but the figure in the generated PDF doesn’t have the dimensions I’d like it to have.
Would anyone be able to provide some insights? Thank you!Hello! I have a slice of a .nii file (essentially a 3D array) I’d like to print out. However, the dimensions of a voxel (a 3D pixel, so basically a cube) aren’t exactly 1mm x 1mm x 1mm. They’re around 1.2mm x 1.1mm x 1.1mm.
I took a slice of the 3D array, so it’s just a 2D image. So, the dimensions of each pixel are ~1.2mm x 1.1mm. I’d like to resize the image such that one pixel is equal to one element of the array is equal to one mm. In theory, it doesn’t seem too difficult. But I just can’t figure this out.
Here’s my code so far:
clear
clc
close all
% Load NIfTI image
info = niftiinfo(‘index_image.nii’);
brain = niftiread(info);
% Extract a slice from the NIfTI image
slice = brain(:, :, 128);
% Create figure
hFig = figure;
% Display the slice with pixel-to-mm scaling
imagesc(slice);
axis off;
axis equal;
colormap(gray)
% Set figure size so that 1 pixel corresponds to 1mm on the printed figure
set(hFig, ‘Units’, ‘centimeters’, ‘Position’, [0 0 24 17.6]);
movegui(hFig, ‘center’);
% Set properties to control the output size
set(hFig, ‘PaperUnits’, ‘centimeters’);
set(hFig, ‘PaperPosition’, [0 0 24 17.6]);
set(hFig, ‘PaperOrientation’, ‘landscape’);
% Export to PDF and open file
print(hFig, ‘-dpdf’, ‘-r0’, ‘out.pdf’);
open(‘out.pdf’);
The reason why use 24 and 17.6 is because my .nii file is 176 x 240 x 256 elements. I’d honestly like the code to work regardless of the size of the array, but even hard coding it into MATLAB isn’t working. The code is compiling just fine, but the figure in the generated PDF doesn’t have the dimensions I’d like it to have.
Would anyone be able to provide some insights? Thank you! Hello! I have a slice of a .nii file (essentially a 3D array) I’d like to print out. However, the dimensions of a voxel (a 3D pixel, so basically a cube) aren’t exactly 1mm x 1mm x 1mm. They’re around 1.2mm x 1.1mm x 1.1mm.
I took a slice of the 3D array, so it’s just a 2D image. So, the dimensions of each pixel are ~1.2mm x 1.1mm. I’d like to resize the image such that one pixel is equal to one element of the array is equal to one mm. In theory, it doesn’t seem too difficult. But I just can’t figure this out.
Here’s my code so far:
clear
clc
close all
% Load NIfTI image
info = niftiinfo(‘index_image.nii’);
brain = niftiread(info);
% Extract a slice from the NIfTI image
slice = brain(:, :, 128);
% Create figure
hFig = figure;
% Display the slice with pixel-to-mm scaling
imagesc(slice);
axis off;
axis equal;
colormap(gray)
% Set figure size so that 1 pixel corresponds to 1mm on the printed figure
set(hFig, ‘Units’, ‘centimeters’, ‘Position’, [0 0 24 17.6]);
movegui(hFig, ‘center’);
% Set properties to control the output size
set(hFig, ‘PaperUnits’, ‘centimeters’);
set(hFig, ‘PaperPosition’, [0 0 24 17.6]);
set(hFig, ‘PaperOrientation’, ‘landscape’);
% Export to PDF and open file
print(hFig, ‘-dpdf’, ‘-r0’, ‘out.pdf’);
open(‘out.pdf’);
The reason why use 24 and 17.6 is because my .nii file is 176 x 240 x 256 elements. I’d honestly like the code to work regardless of the size of the array, but even hard coding it into MATLAB isn’t working. The code is compiling just fine, but the figure in the generated PDF doesn’t have the dimensions I’d like it to have.
Would anyone be able to provide some insights? Thank you! image processing MATLAB Answers — New Questions
Keep underscore symbol using latex as interpreter
Hello All,
I am trying to have one of the entries for my legend (which I like to have latex interpreter for) be the string to the data set I am analyzing (see image). And thus, I want to keep the underscore symbol as is (instead of making what follows a subscript) but still keep the latex as interpreter. I know that making interpeter ‘none’ it will ignore the subscript, but I want to keep the _ symbol in the legend AND use latex interpreter.
I have not figured out a way to do this (in the attached image, latex doesn’t understand the dataName variable so it prints in a different font from the rest of the legend). I want the fonts to match.
Hopefully this makes senseHello All,
I am trying to have one of the entries for my legend (which I like to have latex interpreter for) be the string to the data set I am analyzing (see image). And thus, I want to keep the underscore symbol as is (instead of making what follows a subscript) but still keep the latex as interpreter. I know that making interpeter ‘none’ it will ignore the subscript, but I want to keep the _ symbol in the legend AND use latex interpreter.
I have not figured out a way to do this (in the attached image, latex doesn’t understand the dataName variable so it prints in a different font from the rest of the legend). I want the fonts to match.
Hopefully this makes sense Hello All,
I am trying to have one of the entries for my legend (which I like to have latex interpreter for) be the string to the data set I am analyzing (see image). And thus, I want to keep the underscore symbol as is (instead of making what follows a subscript) but still keep the latex as interpreter. I know that making interpeter ‘none’ it will ignore the subscript, but I want to keep the _ symbol in the legend AND use latex interpreter.
I have not figured out a way to do this (in the attached image, latex doesn’t understand the dataName variable so it prints in a different font from the rest of the legend). I want the fonts to match.
Hopefully this makes sense legend, latex, interpreter MATLAB Answers — New Questions
Closed loop between Simulink and Autoware
Hello,
I am stuck in a problem for a while now, I have created a closed loop between Autoware and Simulink through ROS2. The idea is to analyse the controller node in Autoware. Thats why I am only launching the controller node Controller Node/Orthisnode. This node requires 4 inputs (Trajectory, current position, current acceleration, current steering angle) and one output (controller command). I am publishing these inputs from simulink and subscribing to the output also in simulink.
The controller in Autoware is publishing the output at 32Hz frequency. The problem I am facing is that, controller is not being able to follow the trajectory I am giving from simulink. My guess is the time synchronization problem between Autoware and Simulink. However, I have tried Simulation Pacing option, changed the sample rate of the topics in simulink, also change the fixed step size in Simulink and nothing has been worked.
I am attaching my simulink model and also the ros2 node files where I have changed the message type, hope you can reproduce the results. I thank in advance for any guidance.
Kind regards,
NupurHello,
I am stuck in a problem for a while now, I have created a closed loop between Autoware and Simulink through ROS2. The idea is to analyse the controller node in Autoware. Thats why I am only launching the controller node Controller Node/Orthisnode. This node requires 4 inputs (Trajectory, current position, current acceleration, current steering angle) and one output (controller command). I am publishing these inputs from simulink and subscribing to the output also in simulink.
The controller in Autoware is publishing the output at 32Hz frequency. The problem I am facing is that, controller is not being able to follow the trajectory I am giving from simulink. My guess is the time synchronization problem between Autoware and Simulink. However, I have tried Simulation Pacing option, changed the sample rate of the topics in simulink, also change the fixed step size in Simulink and nothing has been worked.
I am attaching my simulink model and also the ros2 node files where I have changed the message type, hope you can reproduce the results. I thank in advance for any guidance.
Kind regards,
Nupur Hello,
I am stuck in a problem for a while now, I have created a closed loop between Autoware and Simulink through ROS2. The idea is to analyse the controller node in Autoware. Thats why I am only launching the controller node Controller Node/Orthisnode. This node requires 4 inputs (Trajectory, current position, current acceleration, current steering angle) and one output (controller command). I am publishing these inputs from simulink and subscribing to the output also in simulink.
The controller in Autoware is publishing the output at 32Hz frequency. The problem I am facing is that, controller is not being able to follow the trajectory I am giving from simulink. My guess is the time synchronization problem between Autoware and Simulink. However, I have tried Simulation Pacing option, changed the sample rate of the topics in simulink, also change the fixed step size in Simulink and nothing has been worked.
I am attaching my simulink model and also the ros2 node files where I have changed the message type, hope you can reproduce the results. I thank in advance for any guidance.
Kind regards,
Nupur simulink, autoware, simulink pacing, sample rate, subscriber, ros2, publisher MATLAB Answers — New Questions
How to use speechClient(“wav2vec2.0”) command
I need to convert speech into text.I used speechClient("wav2vec2.0") command but got the following error
Error using speechClient (line 26)
Expected apiName to match one of these values:
‘Google’, ‘IBM’, ‘Microsoft’
The input, ‘wav2vec2.0’, did not match any of the valid values.
Please share the solution to resolve the error.I need to convert speech into text.I used speechClient("wav2vec2.0") command but got the following error
Error using speechClient (line 26)
Expected apiName to match one of these values:
‘Google’, ‘IBM’, ‘Microsoft’
The input, ‘wav2vec2.0’, did not match any of the valid values.
Please share the solution to resolve the error. I need to convert speech into text.I used speechClient("wav2vec2.0") command but got the following error
Error using speechClient (line 26)
Expected apiName to match one of these values:
‘Google’, ‘IBM’, ‘Microsoft’
The input, ‘wav2vec2.0’, did not match any of the valid values.
Please share the solution to resolve the error. error using speechclient MATLAB Answers — New Questions
can some one explain me how to pass /give argument to add function instead of void value in simulink function
can some one explain me how to pass /give argument to add function instead of void value in simulink function, How to pass value to add function instead of void value?can some one explain me how to pass /give argument to add function instead of void value in simulink function, How to pass value to add function instead of void value? can some one explain me how to pass /give argument to add function instead of void value in simulink function, How to pass value to add function instead of void value? add function MATLAB Answers — New Questions
Can’t read external drives
I receive the following error when trying to read from several external drives.
"ls: /Volumes/Samsung USB/: Operation not permitted"
I receive similar errors when trying to change to a directory on the external drive.
"Unable to change current folder to ‘/Volumes/Samsung USB’."
The drives are visible in Finder and on Terminal and I can read and write to the drives using other applications. I am running macOS Sonoma version 14.4.1.I receive the following error when trying to read from several external drives.
"ls: /Volumes/Samsung USB/: Operation not permitted"
I receive similar errors when trying to change to a directory on the external drive.
"Unable to change current folder to ‘/Volumes/Samsung USB’."
The drives are visible in Finder and on Terminal and I can read and write to the drives using other applications. I am running macOS Sonoma version 14.4.1. I receive the following error when trying to read from several external drives.
"ls: /Volumes/Samsung USB/: Operation not permitted"
I receive similar errors when trying to change to a directory on the external drive.
"Unable to change current folder to ‘/Volumes/Samsung USB’."
The drives are visible in Finder and on Terminal and I can read and write to the drives using other applications. I am running macOS Sonoma version 14.4.1. file MATLAB Answers — New Questions
Fourier Series,Spectral Analysys
I am trying to implement the spectral analysys for the signal represented in Fig.8. Im trying to do this by using the Fourier Series to generate the signal(it doesnt have to be perfect like in Fig.8) however I am very lost at this point and I dont think thats how the spectral analysys should be looking,i kept trying but I am not able to do it.
Thank you very much for your time.
Sorry for the words not being in english
% Parametrii semnalului dreptunghiular
t1 = 0.2; % Momentul de pornire al creșterii
t2 = 0.8; % Momentul de sfârșit al creșterii
A = 1; % Amplitudinea semnalului dreptunghiular
% Perioada semnalului dreptunghiular (este 1 secundă în acest caz)
T = 1;
% Numărul de armonici pentru seria Fourier (poți ajusta acest număr pentru o aproximare mai bună)
N = 10; % Vom folosi 10 termeni (5 perechi de cosinus și sinus)
% Vectorul de timp
t = linspace(0, T, 1000); % 1000 de puncte între 0 și T pentru o reprezentare mai precisă
% Inițializare semnal dreptunghiular aproximativ
f = zeros(size(t));
% Calcul coeficienți Fourier
a0 = 1; % Coeficientul mediei (pentru semnalul dreptunghiular cu amplitudinea 1)
for n = 1:N
freq = 2*n – 1; % Frecvența armonicului (doar frecvențele impare)
an = (2*A/T) * (1/(pi*freq)) * (sin(2*pi*freq*t2/T) – sin(2*pi*freq*t1/T)); % Coeficientul an
bn = 0; % Pentru semnalul dreptunghiular, coeficientul bn este 0 (nu există termeni sinusoidal)
% Adăugare termen în sumă
f = f + an * cos(2*pi*freq*t);
end
% Adăugare termenul de medie
f = f + a0/2;
% Plotare semnal dreptunghiular aproximativ
figure;
plot(t, f, ‘b’, ‘LineWidth’, 2);
xlabel(‘Timp [secunde]’);
ylabel(‘Amplitudine’);
title(‘Semnal Dreptunghiular Aproximat prin Seria Fourier’);
grid on;
%Analiza spectrala
Y=fft(f);
Pyy=Y.*conj(Y);
s=length(Pyy);
f=20000*(0:s/2-1)/s;
figure;
plot(f,abs(Pyy(1:s/2)));
xlabel(‘Frecventa [Hz]’);
ylabel(‘Spectru’);I am trying to implement the spectral analysys for the signal represented in Fig.8. Im trying to do this by using the Fourier Series to generate the signal(it doesnt have to be perfect like in Fig.8) however I am very lost at this point and I dont think thats how the spectral analysys should be looking,i kept trying but I am not able to do it.
Thank you very much for your time.
Sorry for the words not being in english
% Parametrii semnalului dreptunghiular
t1 = 0.2; % Momentul de pornire al creșterii
t2 = 0.8; % Momentul de sfârșit al creșterii
A = 1; % Amplitudinea semnalului dreptunghiular
% Perioada semnalului dreptunghiular (este 1 secundă în acest caz)
T = 1;
% Numărul de armonici pentru seria Fourier (poți ajusta acest număr pentru o aproximare mai bună)
N = 10; % Vom folosi 10 termeni (5 perechi de cosinus și sinus)
% Vectorul de timp
t = linspace(0, T, 1000); % 1000 de puncte între 0 și T pentru o reprezentare mai precisă
% Inițializare semnal dreptunghiular aproximativ
f = zeros(size(t));
% Calcul coeficienți Fourier
a0 = 1; % Coeficientul mediei (pentru semnalul dreptunghiular cu amplitudinea 1)
for n = 1:N
freq = 2*n – 1; % Frecvența armonicului (doar frecvențele impare)
an = (2*A/T) * (1/(pi*freq)) * (sin(2*pi*freq*t2/T) – sin(2*pi*freq*t1/T)); % Coeficientul an
bn = 0; % Pentru semnalul dreptunghiular, coeficientul bn este 0 (nu există termeni sinusoidal)
% Adăugare termen în sumă
f = f + an * cos(2*pi*freq*t);
end
% Adăugare termenul de medie
f = f + a0/2;
% Plotare semnal dreptunghiular aproximativ
figure;
plot(t, f, ‘b’, ‘LineWidth’, 2);
xlabel(‘Timp [secunde]’);
ylabel(‘Amplitudine’);
title(‘Semnal Dreptunghiular Aproximat prin Seria Fourier’);
grid on;
%Analiza spectrala
Y=fft(f);
Pyy=Y.*conj(Y);
s=length(Pyy);
f=20000*(0:s/2-1)/s;
figure;
plot(f,abs(Pyy(1:s/2)));
xlabel(‘Frecventa [Hz]’);
ylabel(‘Spectru’); I am trying to implement the spectral analysys for the signal represented in Fig.8. Im trying to do this by using the Fourier Series to generate the signal(it doesnt have to be perfect like in Fig.8) however I am very lost at this point and I dont think thats how the spectral analysys should be looking,i kept trying but I am not able to do it.
Thank you very much for your time.
Sorry for the words not being in english
% Parametrii semnalului dreptunghiular
t1 = 0.2; % Momentul de pornire al creșterii
t2 = 0.8; % Momentul de sfârșit al creșterii
A = 1; % Amplitudinea semnalului dreptunghiular
% Perioada semnalului dreptunghiular (este 1 secundă în acest caz)
T = 1;
% Numărul de armonici pentru seria Fourier (poți ajusta acest număr pentru o aproximare mai bună)
N = 10; % Vom folosi 10 termeni (5 perechi de cosinus și sinus)
% Vectorul de timp
t = linspace(0, T, 1000); % 1000 de puncte între 0 și T pentru o reprezentare mai precisă
% Inițializare semnal dreptunghiular aproximativ
f = zeros(size(t));
% Calcul coeficienți Fourier
a0 = 1; % Coeficientul mediei (pentru semnalul dreptunghiular cu amplitudinea 1)
for n = 1:N
freq = 2*n – 1; % Frecvența armonicului (doar frecvențele impare)
an = (2*A/T) * (1/(pi*freq)) * (sin(2*pi*freq*t2/T) – sin(2*pi*freq*t1/T)); % Coeficientul an
bn = 0; % Pentru semnalul dreptunghiular, coeficientul bn este 0 (nu există termeni sinusoidal)
% Adăugare termen în sumă
f = f + an * cos(2*pi*freq*t);
end
% Adăugare termenul de medie
f = f + a0/2;
% Plotare semnal dreptunghiular aproximativ
figure;
plot(t, f, ‘b’, ‘LineWidth’, 2);
xlabel(‘Timp [secunde]’);
ylabel(‘Amplitudine’);
title(‘Semnal Dreptunghiular Aproximat prin Seria Fourier’);
grid on;
%Analiza spectrala
Y=fft(f);
Pyy=Y.*conj(Y);
s=length(Pyy);
f=20000*(0:s/2-1)/s;
figure;
plot(f,abs(Pyy(1:s/2)));
xlabel(‘Frecventa [Hz]’);
ylabel(‘Spectru’); matlab MATLAB Answers — New Questions
Name of Component in Simulink (out.e, out.r)
Does anyone know the component name for out.e, out.r, out.u? I can’t find the component in simulink library.
Thank you.Does anyone know the component name for out.e, out.r, out.u? I can’t find the component in simulink library.
Thank you. Does anyone know the component name for out.e, out.r, out.u? I can’t find the component in simulink library.
Thank you. simulink MATLAB Answers — New Questions
shifting axis on both ends so graph does not have points on extreme end
I have the following attached graph and I need to shift the axis so that I do not have my starting point on the y axis and the last point on the extream end.I have the following attached graph and I need to shift the axis so that I do not have my starting point on the y axis and the last point on the extream end. I have the following attached graph and I need to shift the axis so that I do not have my starting point on the y axis and the last point on the extream end. shifting axis away MATLAB Answers — New Questions
calculate the normal of a 3D plane
Hi! I have a circular plane, whose coordinates of the points of the circumference are:
P = importdata("node_plane.mat");
plot3(P(:,1),P(:,2),P(:,3))
I would like to know if it is correct to calculate the normal of this plane in the following way:
N = cross(P(1,:) – P(2,:), P(3,:) – P(2,:));
N = N/norm(N);Hi! I have a circular plane, whose coordinates of the points of the circumference are:
P = importdata("node_plane.mat");
plot3(P(:,1),P(:,2),P(:,3))
I would like to know if it is correct to calculate the normal of this plane in the following way:
N = cross(P(1,:) – P(2,:), P(3,:) – P(2,:));
N = N/norm(N); Hi! I have a circular plane, whose coordinates of the points of the circumference are:
P = importdata("node_plane.mat");
plot3(P(:,1),P(:,2),P(:,3))
I would like to know if it is correct to calculate the normal of this plane in the following way:
N = cross(P(1,:) – P(2,:), P(3,:) – P(2,:));
N = N/norm(N); normal, plane, 3d MATLAB Answers — New Questions
signature feature extraction matlab code help me
plz helpplz help plz help signature feature extraction matlab code MATLAB Answers — New Questions
Detection of storms from precipitation data
I have 5 min intervals of precipitation data that I am trying to sum into separate storm events.
6 hours of no precipitation (or 72 zero values) signify the end of the storm event.
I’m trying to create a code that
a) Give me the array of matrices (diifferent names e.g M1,M2 etc) of the storm events which lies between 6 hours of no precipitaion.
b) Sums the precip values until the 6 hours of no rainfall and then starts summing the following storm event. The output will be an array with individual storm event depths of the time series.
Kindly help.I have 5 min intervals of precipitation data that I am trying to sum into separate storm events.
6 hours of no precipitation (or 72 zero values) signify the end of the storm event.
I’m trying to create a code that
a) Give me the array of matrices (diifferent names e.g M1,M2 etc) of the storm events which lies between 6 hours of no precipitaion.
b) Sums the precip values until the 6 hours of no rainfall and then starts summing the following storm event. The output will be an array with individual storm event depths of the time series.
Kindly help. I have 5 min intervals of precipitation data that I am trying to sum into separate storm events.
6 hours of no precipitation (or 72 zero values) signify the end of the storm event.
I’m trying to create a code that
a) Give me the array of matrices (diifferent names e.g M1,M2 etc) of the storm events which lies between 6 hours of no precipitaion.
b) Sums the precip values until the 6 hours of no rainfall and then starts summing the following storm event. The output will be an array with individual storm event depths of the time series.
Kindly help. precipitation, storm, matrix, matrix array MATLAB Answers — New Questions