# System of differential equations and fitting with lsqcurvefit

19 ビュー (過去 30 日間)
Monirul Hasan 2019 年 3 月 8 日
コメント済み: darova 2019 年 3 月 18 日
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
%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

#### 2 件のコメント

Alan Weiss 2019 年 3 月 8 日
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 2019 年 3 月 11 日
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.

サインインしてコメントする。

### 回答 (1 件)

Torsten 2019 年 3 月 11 日
function main
%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

#### 2 件のコメント

Monirul Hasan 2019 年 3 月 18 日
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 2019 年 3 月 18 日
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
%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

サインインしてコメントする。

Translated by