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 = 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 

     1

     2

     3

     4

     5

     6

     7

     8

     9

    10

    11

    12

    13

    14

    15

    16

    17

    18

    19

    20

    21

    22

    23

    24

    25

    26

    27

    28

    29

    30

    31

    32

    33

    34

    35

    36

    37

    38

    39

    40

    41

    42

    43

    44

    45

    46

    47

    48

    49

    50

%% 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

Elapsed time is 81.058211 seconds.

%% 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 

[Torsten]

Could you clarify what you think is "not correct" in the results you get from the code above ?

[Katey Stevenson]

Hi @Torsten it's the shape of the solution (i.e., the N_obvs vs laser frequency plot). It is supposed to look like a Lorentzian and instead gives a sort of exponential shape.

Answer 1

[Karan Singh]

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 

     1

     2

     3

     4

     5

     6

     7

     8

     9

    10

    11

    12

    13

    14

    15

    16

    17

    18

    19

    20

    21

    22

    23

    24

    25

    26

    27

    28

    29

    30

    31

    32

    33

    34

    35

    36

    37

    38

    39

    40

    41

    42

    43

    44

    45

    46

    47

    48

    49

    50

%% 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

Elapsed time is 82.378835 seconds.

%% 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