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