Incorrect answer using ode15s (and others) solving ode with time varying coefficients (rate equations); issue with interpolation?
2 views (last 30 days)
Show older comments
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 = 50; %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)
figure(1)
plot(time_interp,A(:,1:p))
figure(2)
plot(time_interp,B(:,1:p))
figure(3)
plot(time_interp,C(:,1:p))
figure(4)
plot(time_interp,D(:,1: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
% 12x1 initial condition array, and a 100x1 time array, there will be
% 100x12 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
2 Comments
Torsten
on 16 May 2024
Could you clarify what you think is "not correct" in the results you get from the code above ?
Answers (1)
Karan Singh
on 18 Jul 2024
Hi Katey,
It seems the main issue lies in the interpolation step within your "makeDiff" function. The interpolation method you are using ('v5cubic') might be causing the coefficients to be interpolated incorrectly, which can affect the ODE solution. Let's try using a different interpolation method such as 'linear', which is often more stable for such problems.
Additionally, ensure that the time points in "time_interp" cover the full range of tspan adequately. Sometimes, inadequate time points can cause interpolation issues.
Here's a slightly modified version of your "solveLIFRateEqns" function with the interpolation method changed to 'linear':
%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 = 50; %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)
figure(1)
plot(time_interp,A(:,1:p))
figure(2)
plot(time_interp,B(:,1:p))
figure(3)
plot(time_interp,C(:,1:p))
figure(4)
plot(time_interp,D(:,1: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
% 12x1 initial condition array, and a 100x1 time array, there will be
% 100x12 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,'linear');
bb = interp1(time_interp,B_col,t,'linear');
cc = interp1(time_interp,C_col,t,'linear');
dd = interp1(time_interp,D_col,t,'linear');
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
Check if this matches your requirement.
Thanks,
Karan
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!