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

조회 수: 3(최근 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
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개)

Community Treasure Hunt

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

Start Hunting!

Translated by