System of differential equations and fitting with lsqcurvefit

Hi, I have a system of equation like this

dN/dt = J(t)/(e*d) - gamma*N.^2, .............. (1) here J(t) is time dependent term- got experimentally

dS/dt = 0.25*gamma*N.^2 - k_rs*s - k_ISC*s + k_RISC*T - Unknown*S.*T - k_SS*S.^2.............(2)

dT/dt = 0.75.*gamma*N.^2 - k_RISC*T+k_ISC*S - k_nT*T................(3)

I have solved eqn(1) first then tried to put the solution [t,y] inside the user defined function BrSolve() by making them Global. It didn't worked.

Please let me know your suggestions t osolve this problem. Feel free to ask me any more details about the problem as well. The Xls data file is also added in the attachment.

clear all, close all,  clc
% loading Exp data
CurrentDensity = xlsread('Modified_Current_density_all','Sheet1','A1:J1001');
Brightness = xlsread('Modified_Brightness_All','Sheet1','A1:J1001');
%declaring solve of eqn (1)/polaron function as global
global t y
%Ydata and xdata for fit
B_70 = Brightness (:,7);
tspan = CurrentDensity (:,1);
% interploating J(t)
Jt = CurrentDensity(:,1);
J = CurrentDensity (:,7);
%solving for polaron number
IC_n = 0;
[t,y] = ode45(@(t,y)Polaron(t,y,J,Jt),tspan,IC_n);
% global t y
figure(3)
plot(t,y)

% Intializing lsqcurvefit
An0 = [1.5e-22]; %sequence [k_rs,k_ISC,k_RISC,k_NRT,k_SS,k_ST,k_TT] ]
lb = [1e-22];% lower bound estimation
ub = [10e-22];
options = optimset('Algorithm', 'trust-region-reflective');   
[An,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@BrSolve,An0,tspan,B_70,lb,ub,options);%sending data for fitting
F = BrSolve(An,tspan);% fitted data saving
% Solving eqn (1)
function n1 = Polaron(t,y,J,Jt)
er=2.9; %relative permitivitty
mobility=1e-5;
d = 15e-7 ; 
e = 1.6e-19*1e3; 
gamma = (e*mobility)/(er*8.8e-12*1e-2);
J = interp1(Jt,J,t);
n1 = J./(e*d)-gamma*y.^2;
end

% function for fit
function Br = BrSolve(An,tspan)
IC_ns = 0;%a
IC_nt = 0;
IC = [IC_ns IC_nt];
[t2,n] = ode45(@(t2,n)Mat(t2,n),tspan,IC);
Br = n(:,1);
%solving eq (2) and eqn (3)
  function dotn = Mat(t2,n)
      k_rs = 1.1000e-02*1e9; 
     k_ISC = 8.7545e-03*1e9; 
     k_RISC = 9.9000e-04*1e9; 
     k_nT = 1.9998e-04*1e9; 
     k_SS = 6.63171736433833e-10*1e9;
     k_ST = 6.59817973391739e-40*1e9;
     k_TT = 7.19219770167270e-24*1e9;
     % gamma part
     er=2.9; %relative permitivitty
     mobility=1e-5;
     d = 15e-7; 
     e = 1.6e-19*1e3; 
     gamma = (e*mobility)/(er*8.8e-12*1e-2);
     y = interp1(tspan,y,t2);
 dotn = zeros(2,1);
     dotn(1) = 0.25*gamma*y.^2-k_rs*n(1)-k_ISC*n(1)+k_RISC*n(2)-An(1)*n(1).*n(2)-k_SS*n(1).^2;%+0.25*k_TT*n(2).^2;
     dotn(2) = 0.75.*gamma*y.^2-k_RISC*n(2)+k_ISC*n(1)-k_nT*n(2);%-1.25*k_TT*n(2).^2;
  end
end

[Alan Weiss]

I tried running your code but got stuck in function dotn = Mat(t2,n). At the line

y = interp1(tspan,y,t2);

I got an error: Undefined function or variable 'y'. Can you please send a correction?

Alan Weiss

MATLAB mathematical toolbox documentation

[Monirul Hasan]

Hi Alan, I believe that's the issue I am having. The [t,y] is the solution of equation 1, which I wanted to get inside the function Mat(t2,n). I tried to make [t, y] as global it didn't worked. I tried to do BrSolve(An,tspan,t,y) it didn't worked. I am not quite sure how to put this inside Mat(t2,n). Please let me know, your suggestion. Thanks.

Answer 1

[Torsten]

function main
  % loading Exp data
  CurrentDensity = xlsread('Modified_Current_density_all','Sheet1','A1:J1001');
  Brightness = xlsread('Modified_Brightness_All','Sheet1','A1:J1001');
  %declaring solve of eqn (1)/polaron function as global
  %Ydata and xdata for fit
  B_70 = Brightness (:,7);
  tspan = CurrentDensity (:,1);
  % interploating J(t)
  Jt = CurrentDensity(:,1);
  J = CurrentDensity (:,7);
  %solving for polaron number
  IC_n = 0;
  [t,y] = ode45(@(t,y)Polaron(t,y,J,Jt),tspan,IC_n);
  figure(3)
  plot(t,y)
  interpfun = @(t2) interp1(t,y,t2);

  % Intializing lsqcurvefit
  An0 = [1.5e-22]; %sequence [k_rs,k_ISC,k_RISC,k_NRT,k_SS,k_ST,k_TT] ]
  lb = [1e-22];% lower bound estimation
  ub = [10e-22];
  options = optimset('Algorithm', 'trust-region-reflective');   
  [An,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@(An,tspan)BrSolve(An,tspan,interpfun),An0,tspan,B_70,lb,ub,options);%sending data for fitting
  F = BrSolve(An,tspan,interpfun);% fitted data saving
end

% Solving eqn (1)
function n1 = Polaron(t,y,J,Jt)
  er=2.9; %relative permitivitty
  mobility=1e-5;
  d = 15e-7 ; 
  e = 1.6e-19*1e3; 
  gamma = (e*mobility)/(er*8.8e-12*1e-2);
  J = interp1(Jt,J,t);
  n1 = J./(e*d)-gamma*y.^2;
end

% function for fit
function Br = BrSolve(An,tspan,interpfun)
  IC_ns = 0;%
  IC_nt = 0;
  IC = [IC_ns IC_nt];
  [t2,n] = ode45(@(t2,n)Mat(t2,n,interpfun),tspan,IC);
  Br = n(:,1);
end
%solving eq (2) and eqn (3)

function dotn = Mat(t2,n,interpfun)
  k_rs = 1.1000e-02*1e9; 
  k_ISC = 8.7545e-03*1e9; 
  k_RISC = 9.9000e-04*1e9; 
  k_nT = 1.9998e-04*1e9; 
  k_SS = 6.63171736433833e-10*1e9;
  k_ST = 6.59817973391739e-40*1e9;
  k_TT = 7.19219770167270e-24*1e9;
  % gamma part
  er=2.9; %relative permitivitty
  mobility=1e-5;
  d = 15e-7; 
  e = 1.6e-19*1e3; 
  gamma = (e*mobility)/(er*8.8e-12*1e-2);
  y = interpfun(t2)

  dotn = zeros(2,1);
  dotn(1) = 0.25*gamma*y.^2-k_rs*n(1)-k_ISC*n(1)+k_RISC*n(2)-An(1)*n(1).*n(2)-k_SS*n(1).^2;%+0.25*k_TT*n(2).^2;
  dotn(2) = 0.75.*gamma*y.^2-k_RISC*n(2)+k_ISC*n(1)-k_nT*n(2);%-1.25*k_TT*n(2).^2;
end

[Monirul Hasan]

Dear Torsten, sorry for my late response. I have tried your modified code. It has some error, I couldn't run. I think 'interpfun' function wasn't defined. Anyway I wil still like your idea of separating Mat() from BrSolve() function. Thanks.

[darova]

Did you try simple variant like this? Take a look at data you have and ode function you try to fit it. Hard to find a solution because of large number. Did you try to pick UKNOWN value manually? Are you sure about your zeros initial conditions?

function main
clc, clear
% loading Exp data
CurrentDensity = xlsread('Modified_Current_density_all');
Brightness = xlsread('Modified_Brightness_All');
%Ydata and xdata for fit
B_70 = Brightness(:,7);
tspan = CurrentDensity(:,1);
J = CurrentDensity(:,7);

An = 1.5e-22;
global iterations
iterations = 0;
% options = odeset('RelTol',1e-4*1e3,...
%                  'AbsTol',[1e-4 1e-4 1e-5]*1e2);
[t,y] = ode45(@(t,y)func(t,y,tspan,J,An),[0 7e-9],[0 0 0]);%, options);

iterations

S = y(:,2);
subplot(211)
plot(t,S,'r')
title('ode')

subplot(212)
plot(tspan,B_70,'b')
title('data')
end



function dy = func(t,y,Jt,Jy,An)
global iterations
k_rs = 1.1000e-02*1e9; 
k_ISC = 8.7545e-03*1e9; 
k_RISC = 9.9000e-04*1e9; 
k_nT = 1.9998e-04*1e9; 
k_SS = 6.63171736433833e-10*1e9;
k_ST = 6.59817973391739e-40*1e9;
k_TT = 7.19219770167270e-24*1e9;
% gamma part
er = 2.9; %relative permitivitty
mobility = 1e-5;
d = 15e-7; 
e = 1.6e-19*1e3; 
gamma = (e*mobility)/(er*8.8e-14);
J = interp1(Jt,Jy,t);

N = y(1);
S = y(2);
T = y(3);
iterations = iterations + 1;
dy = zeros(3,1);
dy(1) = J./(e*d)-gamma*N^2;
dy(2) = 0.25*gamma*N^2 + k_RISC*T - k_ISC*S - k_rs*S - An*S*T - k_SS*S^2;% + 0.25*k_TT*T^2;
dy(3) = 0.75*gamma*N^2 - k_RISC*T + k_ISC*S - k_nT*T;% - 1.25*k_TT*T^2;

end