How to fit data to a PDE using lsqcurvefit/lsqnonlin where y_pdepe is calculated from pdepe output

조회 수: 2 (최근 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).
}
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
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.)
.
Hashim
Hashim 2021년 11월 13일
편집: Hashim 2021년 11월 14일
I guess what I want to do here is to get multiple values of i_num against s_0 and then use lsqcurvefit to estimate k_m. Anyway you think that might be possible please let me know. I have been able to get multiple values of i_num by passing s_0 using arrayfun but do not know how to integrate that with lsqcurvefit so that it automatically gives me an estimate for k_m.
k_m = arrayfun(@PSw_EK_v8, s_bulk, 'un', 0)
Also, I have gone through pretty much all the threads you've mentioned.
For the second part of your suggestion do you mean to say that the equations provided can be solved analytically using symbolic toolbox? because my assumption so far was that the symbolic toolbox is likely to give an significantly long solution (as per an experienced mathematica user). If yes, could you point me towards some resources?

댓글을 달려면 로그인하십시오.

답변 (0개)

카테고리

Help CenterFile Exchange에서 Partial Differential Equation Toolbox에 대해 자세히 알아보기

제품

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by