How to fit data to a PDE using lsqcurvefit/lsqnonlin where y_pdepe is calculated from pdepe output
조회 수: 6 (최근 30일)
이전 댓글 표시
Hi All,
I want to fit some experimental data to my system of PDEs and estimate some system paramters. Trouble is my numerical output is not the direct output of my PDE rather calculated from the PDEPE output (i.e. I am caculating the flux from the PDEPE calculated conc. profile and then using that to calculate current). This is because of my experimental data (current flux vs substrate concentration).
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/799099/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/799104/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/799109/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/799114/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/799119/image.png)
Now my experiemntal data is in the form
s_bulk = [0, 5e-06, 10e-06, 20e-06, 50e-06, 100e-06]
I_exp = [0, 2.35e-03, 2.7e-03, 3.125e-03, 3.29e-03, 3.25e-03]
What I want to do is to fit my pdepe output "I_num" to the "I_exp" given above.
My lsqcurvefit call:
% The purpose of this script is to call PSw_EK_v8 using multiple
% substrate concentrations and calculate the respective current values.
% Using these values we can perhaps try to find the k_m value.
% k_m = [4.92e07, 4.92e07];
s_bulk = [0, 5e-06]; % mol/cm^3
I_exp = [0, 2.35e-03]; % mA/cm^2
% I_num = arrayfun(@PSw_EK_v8, k_m, s_bulk, 'un', 0);
% I_num = cell2mat(I_num);
k_m_0 = 4.92e06;
k_m = lsqcurvefit(@PSw_EK_v8, k_m_0, s_bulk, I_exp);
My pdepe function is given as:
function [I_num] = PSw_EK_v8(k_m, s_bulk)
C.n = 1; % No. of Electrons Involved
C.F = 96485.3329; % C/mol
C.R = 8.314; % J/(mol.K) or CV/(mol.K)
C.Temp = 273.15+37; % K
C.D_s = 7.00E-06; % cm^2/s
C.D_m = 2.81E-06; % cm^2/s (De of polymer)
C.D_e = 0; % because covalently bound
% C.k_m = 4.92E+07; % (mol/cm3)^{-1}s^{-1}
C.kcat = 800; % Turnover No. (1/s)
C.k_e = C.kcat; % Refer to Schnell
C.KM = 2.00E-05; % Michalis Constant (mol/cm^3)
C.P_s = 1; % Partition Coefficient
C.P_m = 1; % Partition Coefficient
% Initial Concentrations
C.e_layer = 5.00E-08; % mol/cm^3
C.m_layer = 1.3E-06; % mol/cm^3
C.Area = 0.0707; % cm^2
C.l = 0.001; % cm
m = 0; % Cartesian Co-ordinates
N = 30; % No. of Points
C.tspan1 = linspace(0, 500, N); % s
% C.tspan1 = logspace(log10(0.000006), log10(500), N); C.tspan1(1) = 0;
xmesh = linspace(0, C.l, N); % cm
% xmesh = logspace(log10(0.000006), log10(C.l), N); xmesh(1) = 0;
% FIRST LINEAR POTENTIAL SWEEP PDEPE Solver Command
E = 0.0; % V
C.E0 = 0.2; % V
C.ScanRt = 0.001; % V/s
C.E_1 = E+(C.ScanRt*C.tspan1); % Potential Sweep 1
C.epsilon1 = ((C.n*C.F)/(C.R*C.Temp)).*(C.E_1-C.E0);
C.c_mred1 = C.m_layer./(1+exp(C.epsilon1)); % Mred
% FIRST POTENTIAL SWEEP Initial & Boundary Condition Call
ic_arg = {@(x)s_bulk*ones(size(N)) ; @(x)C.e_layer*ones(size(N)); ...
@(x)C.m_layer*ones(size(N))};
IC = @(x)PDE_PSw_EK_IC(x, ic_arg);
BC = @(xl, cl, xr, cr, t)PDE_PSw_EK_BC(xl, cl, xr, cr, t, C, s_bulk);
% Calling pdepe solver after time mesh refinement
% MaxStep sets an upper bound on the size of any step taken by the solver.
% If the equation has periodic behavior, for example, then setting MaxStep
% to a fraction of the period ensures that the solver does not enlarge the
% step so much that it steps over an area of interest.
optns = odeset('MaxStep',1e-01,'RelTol',1e-09,'AbsTol',1e-10);
sol1 = pdepe(m, @(x, t, c, DcDx)PDE_PSw_EK(x, t, c, DcDx, C, k_m), ...
IC, BC, xmesh, C.tspan1, optns);
% Concentration Profiles c(i, j, k)(Solutions)
c1 = sol1(:, :, 1); % Substrate Conc.
c2 = sol1(:, :, 2); % Ered Conc.
c3 = sol1(:, :, 3); % Mred Conc.
c4 = C.m_layer-c3; % Mox Conc.
% Plot Species Conc. in the Layer for FIRST POTENTIAL SWEEP
figure(1);
[hAx,hLine1,hLine2]= plotyy(xmesh,c4(N,:),xmesh,c2(N,:));
xlabel('x / cm');
ylabel(hAx(1),'[c_{m_{ox}}] / mol cm^{-^3}') % left y-axis
ylabel(hAx(2),'[c_{e_{red}}] / mol cm^{-^3}') % right y-axis
title('Concentratrions for FIRST POTENTIAL SWEEP');
% % Plot Species Conc. in the Layer
% figure(11)
% plot(xmesh, c1(N,:));
% xlabel('x / cm'); ylabel('substrate [s]');
% [~, dc1Sdx_0(counter)] = pdeval(m, xmesh, c1(counter,:), 0); % Substrate Flux at x=0
% [~, dc1Sdx_l(counter)] = pdeval(m, xmesh, c1(counter,:),C.l); % Substrate Flux at x=d
[~, dc4Mdx_0] = pdeval(m, xmesh, c4(N,:), 0); % Mox Flux at x=0
% [~, dc4Mdx_l(counter)] = pdeval(m, xmesh, c4(counter,:),C.l); % Mox Flux at x=d
j_num_ano = (-C.D_m .* dc4Mdx_0);
I_num = (C.n * C.F .* j_num_ano)/(C.Area);
% I vs E for FIRST POTENTIAL SWEEP
figure(4);
plot(C.E_1, I_num, 'r--', 'LineWidth',0.5);
xlabel('E / V'); ylabel('I / A');
hold on;
z = 1;
end
%% pdepe Function PSw_EK_v8
%%
% Notes C.D_e is 0 because enzyme is immobilized in the polymeric matrix
% and cannot diffuse across the layer.
function [cc, ff, ss] = PDE_PSw_EK(x, t, c, DcDx, C, k_m)
% Substrate; Ered; Mred;
cc = [ 1; 1; 1];
ff = [C.D_s; C.D_e; C.D_m].*DcDx;
menten = (C.k_e*c(1))*(C.e_layer-c(2))/(C.KM+c(1));
mred = k_m*(C.m_layer-c(3))*c(2);
ss = [-menten; -mred+menten; +mred];
end
%% Initial Condition --> Initial Concentrations of Species
%%
function c0 = PDE_PSw_EK_IC(x, ic_arg)
% Substrate; Ered; Mred;
c0 = [ic_arg{1}(x); ic_arg{2}(x); ic_arg{3}(x)];
end
%% Boundary Condition - 1
%%
function [pl, ql, pr, qr] = PDE_PSw_EK_BC(xl, cl, xr, cr, t, C, s_bulk)
% ---------------------
% | | |
% | | |
% | | |
% | | |
% | | |
% | | |
% | | |
% -------|-------------
% pl pr
% Order --> Substrate; Mox; Mred; Eox; Ered
% Substrate;Er; Mox;
pl = [0 ; 0; cl(3)-interp1(C.tspan1,C.c_mred1,t)];
ql = [1 ; 1; 0];
pr = [cr(1)-s_bulk ; 0; 0];
qr = [0 ; 1; 1];
end
댓글 수: 2
Star Strider
2021년 11월 13일
I never tried to do this with pdepe (and I have very little experienced with pdepe) however a typical approach for ordinary differential equations is presented in Parameter Estimation for a System of Differential Equations and similar threads.
If the partial differential equations are linear and have constant coefficients (they appear to, however I could be missing something), another option is to use the Symbolic Math Toolbox to first take the Laplace transform of the time derivatives, creating a system of ordinary differential equations in the spatial (or other) derivatives that can then be solved analytically. Then invert them back to the time domain and fit the data with that. (That was a grad school final exam problem, and something of a challenge to do my hand, however I got it right.)
.
답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Partial Differential Equation Toolbox에 대해 자세히 알아보기
제품
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!