System of differential equations and fitting with lsqcurvefit
1 回表示 (過去 30 日間)
古いコメントを表示
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
回答 (1 件)
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 件のコメント
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
参考
カテゴリ
Help Center および File Exchange で Linear Algebra についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!