MATLAB Answers

0

System of differential equations and fitting with lsqcurvefit

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

  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.

サインイン to comment.

1 件の回答

Torsten
回答者: Torsten
2019 年 3 月 11 日

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

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

サインイン to comment.



Translated by