How to input a novel boundary condition for a coupled PDE system

Question

Hashim 2024-8-14，22:59

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2145369-how-to-input-a-novel-boundary-condition-for-a-coupled-pde-system

评论： Torsten 2024-8-20，13:56

Hi all,

I am trying to simulate a coupled PDE system with a non typical boundary conditions using the pde1dm solver which is based on the pdepe solver in MATLAB.

My coupled PDE system is as such.

The system is solved for two time periods

For first time period

For s at the L.H.S we have

For s at the R.H.S we have

For species m at the L.H.S we have the to set m value to a constant

For species m at the R.H.S we have the no flux condition

For second time period

For s at the L.H.S and R.H.S we have the same such as

Also for the m R.H.S we have the same boundary condition.

But, for the species

we have the L.H.S boundary condition for timr period

. It is this boundary condtiuon that I am seeking help with. How can I go about implementing this on the L.H.S for

.

function [c25, c2, vode] = pde1dm_PS_OCP_v1(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_r)
% Piecewise PDE Implementation 
% Pulsing Simulation 
% Constants Used in the Nernst Equation   
n           = 2;                                             % No. of Electrons Involved
F           = 96485.3329;                                    % C/mol
R           = 8.314;                                         % J/(mol.K) or CV/(mol.K)
Temp        = 273.15+20;                                     % K
N           = 21; 
m           = 0; 
Area        = 0.0011;                                                     % Geometric Area (cm^2)
m_layer     = 4.13E-05;                                                   % m_layer (mol/cm^3)     
D_m         = 8.14E-08;                                                   % cm^2/s (De of polymer)
k           = n*F*Area*l*m_layer;
% DIMENSIONLESS PARAMETERS
chi         = linspace(0, 1, N);
tau         = linspace(0, tau_M1, N);                               % Dimensionless Time 
% FIRST POTENTIAL STEP    
E0          = 0.22;                                                 
E_set       = .350;                                      
epsilon1    = ((n*F)/(R*Temp))*(E_set-E0);                          % Dimensionless Potential 
c_M_ox      = 1/(1+exp(-epsilon1));                                 % Mox BC 
ic_arg_1    = {@(x)ones(size(N)) ; @(x)zeros(size(N))};
IC          =  @(x)PDE_ND_PS_EK_IC(x, ic_arg_1); 
BC          =  @(xl, yl, xr, yr, t)PDE_PS_EK_BC(xl, yl, xr, yr, c_M_ox); 
sol1        = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r), ... 
                    IC, BC, chi, tau);
 
c1          = sol1(:, :, 1);                                        % Substrate Conc.  
c2          = sol1(:, :, 2);                                        % Mox Conc.
               
% OPEN CIRCUIT POTENTIAL    
tau2         = linspace(tau_M1, tau_M2, N);                         % Dimensionless Time 
ic_arg_1     = {@(x)interp1(chi, sol1(N, :, 1), x) ; @(x)interp1(chi, sol1(N, :, 2), x)};
IC2          =  @(x)PDE_ND_PS_EK_IC(x, ic_arg_1); 
BC2          =  @(xl, yl, xr, yr, t, v)PDE_PSw_EK_BC_2(xl, yl, xr, yr, v); 
xOde         =  [1, .9].';
ode_IC       =  @() ode_ICFun(tau_M1);
opts.vectorized='on';
% sol2         = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(x, t, c, DcDx, kappa, eta, gamma, mu, D_r), ... 
%                     IC2, BC2, chi, tau2);
[sol2, vode] = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(x, t, c, DcDx, kappa, eta, gamma, mu, D_r), ... 
               IC2, BC2, chi, tau2, @odeFunc, ode_IC, xOde(2));
               
% Concentration Profiles c(i, j, k)(Solutions) 
c14          = [c1; sol2(:, :, 1)];                                 % Substrate Conc.  
c25          = [c2; sol2(:, :, 2)];                                 % Mox Conc. 
c36          = 1-c25;                                               % Mred Conc. 
mox           = abs(sol2(:, 1, 2));
mred          = 1-mox;
E_PS          = E_set.*ones(1,N); 
for counter2  = 1:N
    E_OCP(counter2) = E0 + (((R*Temp)/(n*F).*(log(mox(counter2,:)./mred(counter2,:)))));
end
E_array       = [E_PS, E_OCP]; 
figure(2);
plot(chi, c25(N,:));
xlabel('\chi'); ylabel('m_{ox}');
function [cc, ff, ss] = PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r)
 
    %                     S;          Mox;   
    cc              = [ D_r;           1]; 
    
    ff              = [ 1;             1].*DcDx;
    
    S_kin           = ((gamma/eta)*(kappa^2)*c(1)*c(2))./...
                      (gamma.*c(2).*(1+(mu.*c(1)))+c(1));
    M_kin           = ((kappa^2).*c(1).*c(2))./...
                      (gamma.*c(2).*(1+(mu.*c(1)))+c(1));
    ss              = [-S_kin;  -M_kin]; 
           
end
function c0     = PDE_ND_PS_EK_IC(x, ic_arg_1)
   
    %                         Substrate;                Mox;  
    c0              = [ic_arg_1{1}(x).';  ic_arg_1{2}(x).'];
end
% % % % % % % % % % % % % % % % For 0<tau_m<tau_1
function [pl, ql, pr, qr] = PDE_PS_EK_BC(xl, yl, xr, yr, c_M_ox)
% ---------------------
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% -------|-------------
%        pl           pr 
    %                Substrate;         Mediator;
    pl     = [0               ;    yl(2)-c_M_ox];     
    
    ql     = [1               ;               0];
    
    pr     = [yr(1)-1         ;               0];
    
    qr     = [0               ;               1]; 
end
% % % % % % % % % % % % % % % % For tau_1<tau_m<tau_2
function [pl, ql, pr, qr] = PDE_PSw_EK_BC_2(xl, yl, xr, yr, t, v, vdot) 
    pl     = [0               ;               0];     
    
    ql     = [1               ;               1];
    
    pr     = [yr(1)-1         ;               0];
    
    qr     = [0               ;               1]; 
end
function f = odeFunc(x, t, c, DcDx, v, vdot)
    f = (DcDx*(1E-03./(v(1)*(1-v(1)))))-vdot(1);
end
function v0=ode_ICFun(tau_0)
    v0 = ones(tau_0);
end
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-8-14，23:44

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2145369-how-to-input-a-novel-boundary-condition-for-a-coupled-pde-system#answer_1499134

在 MATLAB Online 中打开

xOde = [0 1].';  % Define coupling points for ODE definition
function [pl, ql, pr, qr] = PDE_PSw_EK_BC_2(xl, yl, xr, yr, t, v, vdot) 
    ...
    pl(2) = yl(2) - v;
    ql(2) = 0;
    ...
end
function f = odeFunc(x, t, c, DcDx, v, vdot)
  % Assuming first index in DcDx and c is equation, second index is coupling point in xOde vector 
  % (might be vice versa, look up in the User's Guide of pde1dm)
  f = DcDx(2,1) - k/(c(2,2)*(1-c(2,2))*vdot     
end
function v0=ode_ICFun(tau_0)
  v0 = m_ox %@t=0,chi=0
end

18 个评论
显示 16更早的评论隐藏 16更早的评论

Hashim 2024-8-15，10:41

在 MATLAB Online 中打开

pde1dM_manual.pdf

Hi Torsten, thank you for replying. So, for a bit of backgriound I am using the pde1dm function which I believe is written by Mr. Greene (tutorial attached). This function can accomodate pdes and odes simultaneously. So for the left hand side I have the boundary given as. I am solving for two concentration

and s. It is the

for which I need to apply the novel boundary condition.

But it is not a separate ODE but a boundary condition so the

is in fact the same c as defined in the pdefun function. It has the value=1 at time t=0 i.e. initial condition is 1. But i suspect you need to solve the ODE and then apply it as boundary condition using the pdefun. So far as you can see I have not done that rather as a placeholder I have just applied the no flux condition on the LHS for

. So right now I am trying to solve for "v" in the odeFunc and I do get an array but I am not sure if this is the right approach and then there is the question of trying to apply it in the PDE_PSw_EK_BC_2 function.

function [pl, ql, pr, qr] = PDE_PSw_EK_BC_2(xl, yl, xr, yr, t, v, vdot) 
    pl     = [0               ;               0];     
    
    ql     = [1               ;               1];
    
    pr     = [yr(1)-1         ;               0];
    
    qr     = [0               ;               1]; 
end

Hashim 2024-8-15，13:23

编辑：Hashim 2024-8-15，15:00

在 MATLAB Online 中打开

Hi torsten,

Well, I've been trying and so far it seems the ode results are being applied to the LHS if you look at c25 as it is the entire conc. profile for m_ox. There are problems though.

function [c25, c14, c2, vode] = pde1dm_PS_OCP_v1(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_r)
% Piecewise PDE Implementation 
% Pulsing Simulation 
% Constants Used in the Nernst Equation   
n           = 2;                                             % No. of Electrons Involved
F           = 96485.3329;                                    % C/mol
R           = 8.314;                                         % J/(mol.K) or CV/(mol.K)
Temp        = 273.15+20;                                     % K
N           = 21; 
m           = 0; 
Area        = 0.0011;                                                     % Geometric Area (cm^2)
m_layer     = 4.13E-05;                                                   % m_layer (mol/cm^3)     
D_m         = 8.14E-08;                                                   % cm^2/s (De of polymer)
k           = n*F*Area*l*m_layer;
% DIMENSIONLESS PARAMETERS
chi         = linspace(0, 1, N);
tau         = linspace(0, tau_M1, N);                               % Dimensionless Time 
% FIRST POTENTIAL STEP    
E0          = 0.22;                                                 
E_set       = .350;                                      
epsilon1    = ((n*F)/(R*Temp))*(E_set-E0);                          % Dimensionless Potential 
c_M_ox      = 1/(1+exp(-epsilon1));                                 % Mox BC 
ic_arg_1    = {@(x)ones(size(N)) ; @(x)zeros(size(N))};
IC          =  @(x)PDE_ND_PS_EK_IC(x, ic_arg_1); 
BC          =  @(xl, yl, xr, yr, t)PDE_PS_EK_BC(xl, yl, xr, yr, c_M_ox); 
sol1        = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r), ... 
                    IC, BC, chi, tau);
 
c1          = sol1(:, :, 1);                                        % Substrate Conc.  
c2          = sol1(:, :, 2);                                        % Mox Conc.
               
% OPEN CIRCUIT POTENTIAL    
tau2         = linspace(tau_M1, tau_M2, N);                         % Dimensionless Time 
ic_arg_1     = {@(x)interp1(chi, sol1(N, :, 1), x) ; @(x)interp1(chi, sol1(N, :, 2), x)};
IC2          =  @(x)PDE_ND_PS_EK_IC(x, ic_arg_1); 
BC2          =  @(xl, yl, xr, yr, t, v)PDE_PSw_EK_BC_2(xl, yl, xr, yr, v); 
xOde         =  [0, 1].';
ode_IC       =  @() ode_ICFun(tau_M1);
opts.vectorized='on';
% sol2         = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(x, t, c, DcDx, kappa, eta, gamma, mu, D_r), ... 
%                     IC2, BC2, chi, tau2);
[sol2, vode] = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(x, t, c, DcDx, kappa, eta, gamma, mu, D_r), ... 
               IC2, @PDE_PSw_EK_BC_2, chi, tau2, @odeFunc, ode_IC, .9);
               
% Concentration Profiles c(i, j, k)(Solutions) 
c14          = [c1; sol2(:, :, 1)];                                 % Substrate Conc.  
c25          = [c2; sol2(:, :, 2)];                                 % Mox Conc. 
c36          = 1-c25;                                               % Mred Conc. 
mox           = abs(sol2(:, 1, 2));
mred          = 1-mox;
E_PS          = E_set.*ones(1,N); 
for counter2  = 1:N
    E_OCP(counter2) = E0 + (((R*Temp)/(n*F).*(log(mox(counter2,:)./mred(counter2,:)))));
end
E_array       = [E_PS, E_OCP]; 
figure(2);
plot(chi, c25(N,:));
xlabel('\chi'); ylabel('m_{ox}');
function [cc, ff, ss] = PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r)
 
    %                     S;          Mox;   
    cc              = [ D_r;           1]; 
    
    ff              = [ 1;             1].*DcDx;
    
    S_kin           = ((gamma/eta)*(kappa^2)*c(1)*c(2))./...
                      (gamma.*c(2).*(1+(mu.*c(1)))+c(1));
    M_kin           = ((kappa^2).*c(1).*c(2))./...
                      (gamma.*c(2).*(1+(mu.*c(1)))+c(1));
    ss              = [-S_kin;  -M_kin]; 
           
end
function c0     = PDE_ND_PS_EK_IC(x, ic_arg_1)
   
    %                         Substrate;                Mox;  
    c0              = [ic_arg_1{1}(x).';  ic_arg_1{2}(x).'];
end
% % % % % % % % % % % % % % % % For 0<tau_m<tau_1
function [pl, ql, pr, qr] = PDE_PS_EK_BC(xl, yl, xr, yr, c_M_ox)
% ---------------------
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% -------|-------------
%        pl           pr 
    %                Substrate;         Mediator;
    pl     = [0               ;    yl(2)-c_M_ox];     
    
    ql     = [1               ;               0];
    
    pr     = [yr(1)-1         ;               0];
    
    qr     = [0               ;               1]; 
end
% % % % % % % % % % % % % % % % For tau_1<tau_m<tau_2
function [pl, ql, pr, qr] = PDE_PSw_EK_BC_2(xl, yl, xr, yr, t, v, vdot) 
    pl     = [0               ;               yl(2)-v];     
    
    ql     = [1               ;               0];
    
    pr     = [yr(1)-1         ;               0];
    
    qr     = [0               ;               1]; 
end
function f = odeFunc(x, t, c, DcDx, v, vdot)
    f = (DcDx*(1E03./(v(1)*(1-v(1)))))-vdot(1);
%     f=DcDx(2,1) - 1/(c(2,2)*(1-c(2,2)))*vdot  ;
end
function v0=ode_ICFun(tau_0)
    v0 = ones(tau_0);
end
end

Hashim 2024-8-15，16:06

编辑：Torsten 2024-8-15，19:32

在 MATLAB Online 中打开

So, here is the script where I've made all the changes I believe that you suggested. It gves the following error.

[c25, c14, c2, vode] = pde1dm_PS_OCP_v1(1, 1, 1, 1, 1, 2, 1, 1)

v0 = 1

c25 = 42x21

1.0000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0000 0.8727 0.7492 0.6321 0.5238 0.4260 0.3398 0.2657 0.2036 0.1527 0.1121 0.0806 0.0566 0.0389 0.0261 0.0171 0.0110 0.0070 0.0046 0.0033 0.0028 1.0000 0.9078 0.8175 0.7300 0.6463 0.5670 0.4930 0.4247 0.3623 0.3062 0.2564 0.2127 0.1749 0.1428 0.1160 0.0942 0.0769 0.0639 0.0548 0.0494 0.0477 1.0000 0.9233 0.8481 0.7747 0.7038 0.6359 0.5712 0.5103 0.4535 0.4009 0.3527 0.3092 0.2703 0.2361 0.2066 0.1818 0.1616 0.1460 0.1349 0.1282 0.1260 1.0000 0.9329 0.8671 0.8028 0.7405 0.6805 0.6231 0.5686 0.5173 0.4695 0.4252 0.3848 0.3483 0.3159 0.2876 0.2636 0.2439 0.2286 0.2176 0.2110 0.2088 1.0000 0.9400 0.8812 0.8239 0.7682 0.7145 0.6630 0.6141 0.5679 0.5246 0.4845 0.4477 0.4143 0.3846 0.3587 0.3365 0.3183 0.3041 0.2939 0.2878 0.2857 1.0000 0.9458 0.8928 0.8410 0.7908 0.7424 0.6960 0.6519 0.6102 0.5711 0.5348 0.5015 0.4713 0.4443 0.4207 0.4006 0.3840 0.3711 0.3618 0.3562 0.3543 1.0000 0.9508 0.9026 0.8556 0.8101 0.7663 0.7242 0.6843 0.6465 0.6111 0.5783 0.5481 0.5207 0.4963 0.4749 0.4567 0.4417 0.4299 0.4215 0.4164 0.4147 1.0000 0.9550 0.9110 0.8682 0.8268 0.7869 0.7486 0.7123 0.6780 0.6459 0.6161 0.5887 0.5639 0.5417 0.5223 0.5057 0.4921 0.4814 0.4738 0.4692 0.4676 1.0000 0.9587 0.9183 0.8791 0.8412 0.8048 0.7699 0.7367 0.7054 0.6762 0.6490 0.6241 0.6015 0.5813 0.5637 0.5486 0.5362 0.5265 0.5195 0.5154 0.5140

<mw-icon class=""></mw-icon>

c14 = 42x21

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 0.9893 0.9894 0.9898 0.9905 0.9913 0.9922 0.9932 0.9942 0.9952 0.9961 0.9969 0.9977 0.9983 0.9987 0.9991 0.9994 0.9996 0.9998 0.9999 0.9999 1.0000 0.9788 0.9789 0.9793 0.9800 0.9809 0.9819 0.9831 0.9845 0.9859 0.9873 0.9887 0.9901 0.9915 0.9928 0.9940 0.9951 0.9962 0.9971 0.9981 0.9990 1.0000 0.9683 0.9685 0.9689 0.9696 0.9705 0.9717 0.9730 0.9745 0.9761 0.9778 0.9796 0.9814 0.9833 0.9852 0.9872 0.9891 0.9911 0.9932 0.9953 0.9976 1.0000 0.9580 0.9582 0.9586 0.9594 0.9604 0.9616 0.9630 0.9647 0.9665 0.9685 0.9706 0.9729 0.9753 0.9778 0.9805 0.9832 0.9862 0.9893 0.9926 0.9961 1.0000 0.9481 0.9483 0.9488 0.9495 0.9506 0.9519 0.9535 0.9554 0.9574 0.9597 0.9622 0.9649 0.9678 0.9710 0.9743 0.9779 0.9817 0.9857 0.9901 0.9949 1.0000 0.9387 0.9389 0.9394 0.9402 0.9414 0.9429 0.9446 0.9467 0.9490 0.9516 0.9545 0.9576 0.9611 0.9648 0.9687 0.9730 0.9777 0.9826 0.9880 0.9938 1.0000 0.9298 0.9300 0.9306 0.9315 0.9328 0.9344 0.9364 0.9387 0.9413 0.9442 0.9474 0.9510 0.9549 0.9592 0.9638 0.9687 0.9741 0.9799 0.9861 0.9928 1.0000 0.9216 0.9218 0.9224 0.9234 0.9248 0.9266 0.9288 0.9313 0.9342 0.9374 0.9410 0.9450 0.9494 0.9541 0.9593 0.9649 0.9709 0.9774 0.9844 0.9919 1.0000 0.9141 0.9143 0.9149 0.9160 0.9175 0.9195 0.9218 0.9245 0.9277 0.9312 0.9352 0.9396 0.9444 0.9496 0.9553 0.9615 0.9681 0.9752 0.9829 0.9912 1.0000

<mw-icon class=""></mw-icon>

c2 = 21x21

1.0000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0000 0.8727 0.7492 0.6321 0.5238 0.4260 0.3398 0.2657 0.2036 0.1527 0.1121 0.0806 0.0566 0.0389 0.0261 0.0171 0.0110 0.0070 0.0046 0.0033 0.0028 1.0000 0.9078 0.8175 0.7300 0.6463 0.5670 0.4930 0.4247 0.3623 0.3062 0.2564 0.2127 0.1749 0.1428 0.1160 0.0942 0.0769 0.0639 0.0548 0.0494 0.0477 1.0000 0.9233 0.8481 0.7747 0.7038 0.6359 0.5712 0.5103 0.4535 0.4009 0.3527 0.3092 0.2703 0.2361 0.2066 0.1818 0.1616 0.1460 0.1349 0.1282 0.1260 1.0000 0.9329 0.8671 0.8028 0.7405 0.6805 0.6231 0.5686 0.5173 0.4695 0.4252 0.3848 0.3483 0.3159 0.2876 0.2636 0.2439 0.2286 0.2176 0.2110 0.2088 1.0000 0.9400 0.8812 0.8239 0.7682 0.7145 0.6630 0.6141 0.5679 0.5246 0.4845 0.4477 0.4143 0.3846 0.3587 0.3365 0.3183 0.3041 0.2939 0.2878 0.2857 1.0000 0.9458 0.8928 0.8410 0.7908 0.7424 0.6960 0.6519 0.6102 0.5711 0.5348 0.5015 0.4713 0.4443 0.4207 0.4006 0.3840 0.3711 0.3618 0.3562 0.3543 1.0000 0.9508 0.9026 0.8556 0.8101 0.7663 0.7242 0.6843 0.6465 0.6111 0.5783 0.5481 0.5207 0.4963 0.4749 0.4567 0.4417 0.4299 0.4215 0.4164 0.4147 1.0000 0.9550 0.9110 0.8682 0.8268 0.7869 0.7486 0.7123 0.6780 0.6459 0.6161 0.5887 0.5639 0.5417 0.5223 0.5057 0.4921 0.4814 0.4738 0.4692 0.4676 1.0000 0.9587 0.9183 0.8791 0.8412 0.8048 0.7699 0.7367 0.7054 0.6762 0.6490 0.6241 0.6015 0.5813 0.5637 0.5486 0.5362 0.5265 0.5195 0.5154 0.5140

<mw-icon class=""></mw-icon>

vode = 21x1

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

<mw-icon class=""></mw-icon>

The scipt is as:

function [c25, c14, c2, vode] = pde1dm_PS_OCP_v1(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_r)
% Piecewise PDE Implementation 
% Pulsing Simulation 
% Constants Used in the Nernst Equation   
n           = 2;                                             % No. of Electrons Involved
F           = 96485.3329;                                    % C/mol
R           = 8.314;                                         % J/(mol.K) or CV/(mol.K)
Temp        = 273.15+20;                                     % K
N           = 21; 
m           = 0; 
Area        = 0.0011;                                                     % Geometric Area (cm^2)
m_layer     = 4.13E-05;                                                   % m_layer (mol/cm^3)     
D_m         = 8.14E-08;                                                   % cm^2/s (De of polymer)
k           = n*F*Area*l*m_layer;
% DIMENSIONLESS PARAMETERS
chi         = linspace(0, 1, N);
tau         = linspace(0, tau_M1, N);                               % Dimensionless Time 
% FIRST POTENTIAL STEP    
E0          = 0.22;                                                 
E_set       = .350;                                      
epsilon1    = ((n*F)/(R*Temp))*(E_set-E0);                          % Dimensionless Potential 
c_M_ox      = 1/(1+exp(-epsilon1));                                 % Mox BC 
ic_arg_1    = {@(x)ones(size(N)) ; @(x)zeros(size(N))};
IC          =  @(x)PDE_ND_PS_EK_IC(x, ic_arg_1); 
BC          =  @(xl, yl, xr, yr, t)PDE_PS_EK_BC(xl, yl, xr, yr, c_M_ox); 
sol1        = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r), ... 
                    IC, BC, chi, tau);
 
c1          = sol1(:, :, 1);                                        % Substrate Conc.  
c2          = sol1(:, :, 2);                                        % Mox Conc.
               
% OPEN CIRCUIT POTENTIAL    
tau2         = linspace(tau_M1, tau_M2, N);                         % Dimensionless Time 
ic_arg_1     = {@(x)interp1(chi, sol1(N, :, 1), x) ; @(x)interp1(chi, sol1(N, :, 2), x)};
IC2          =  @(x)PDE_ND_PS_EK_IC(x, ic_arg_1); 
BC2          =  @(xl, yl, xr, yr, t, v)PDE_PSw_EK_BC_2(xl, yl, xr, yr, v); 
xOde         =  0;
ode_IC       =  @() ode_ICFun(tau_M1);
opts.vectorized='on';
% sol2         = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(x, t, c, DcDx, kappa, eta, gamma, mu, D_r), ... 
%                     IC2, BC2, chi, tau2);
[sol2, vode] = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(x, t, c, DcDx, kappa, eta, gamma, mu, D_r), ... 
               IC2, @PDE_PSw_EK_BC_2, chi, tau2, @odeFunc, ode_IC, xOde);
               
% Concentration Profiles c(i, j, k)(Solutions) 
c14          = [c1; sol2(:, :, 1)];                                 % Substrate Conc.  
c25          = [c2; sol2(:, :, 2)];                                 % Mox Conc. 
c36          = 1-c25;                                               % Mred Conc. 
mox           = abs(sol2(:, 1, 2));
mred          = 1-mox;
E_PS          = E_set.*ones(1,N); 
for counter2  = 1:N
    E_OCP(counter2) = E0 + (((R*Temp)/(n*F).*(log(mox(counter2,:)./mred(counter2,:)))));
end
E_array       = [E_PS, E_OCP]; 
% This is the case finder function
% Video(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_m, N, chi, c14, c25, c36, E_array)
figure(2);
plot(chi, c25(end,:));
xlabel('\chi'); ylabel('m_{ox}');
function [cc, ff, ss] = PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r)
 
    %                     S;          Mox;   
    cc              = [ D_r;           1]; 
    
    ff              = [ 1;             1].*DcDx;
    
    S_kin           = ((gamma/eta)*(kappa^2)*c(1)*c(2))./...
                      (gamma.*c(2).*(1+(mu.*c(1)))+c(1));
    M_kin           = ((kappa^2).*c(1).*c(2))./...
                      (gamma.*c(2).*(1+(mu.*c(1)))+c(1));
    ss              = [-S_kin;  -M_kin]; 
           
end
function c0     = PDE_ND_PS_EK_IC(x, ic_arg_1)
   
    %                         Substrate;                Mox;  
    c0              = [ic_arg_1{1}(x).';  ic_arg_1{2}(x).'];
end
% % % % % % % % % % % % % % % % For 0<tau_m<tau_1
function [pl, ql, pr, qr] = PDE_PS_EK_BC(xl, yl, xr, yr, c_M_ox)
% ---------------------
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% -------|-------------
%        pl           pr 
    %                Substrate;         Mediator;
    pl     = [0               ;    yl(2)-c_M_ox];     
    
    ql     = [1               ;               0];
    
    pr     = [yr(1)-1         ;               0];
    
    qr     = [0               ;               1]; 
end
% % % % % % % % % % % % % % % % For tau_1<tau_m<tau_2
function [pl, ql, pr, qr] = PDE_PSw_EK_BC_2(xl, yl, xr, yr, t, v, vdot) 
    pl     = [0               ;               yl(2)-v];     
    
    ql     = [1               ;               0];
    
    pr     = [yr(1)-1         ;               0];
    
    qr     = [0               ;               1]; 
end
function f = odeFunc(t, v, vdot, x, c, DcDx)
%     f = (DcDx*(1E03./(v(1)*(1-v(1)))))-vdot(1);
    f=DcDx(2) - 10/(c(2)*(1-c(2)))*vdot(1);
end
function v0=ode_ICFun(tau_0)
    v0 = ones(tau_0)
end
end

Hashim 2024-8-15，20:21

编辑：Hashim 2024-8-15，20:26

在 MATLAB Online 中打开

Hi, here is the code I am running. For me I can see the v array changing and I can also see it is applied as a boundary if you see the c25. Can you see any mistakes in this script?

BTW Thank you so very much guys!

[c25, c2, vode, chi] = pde1dm_PS_OCP_v1(1, 1, 0.2, 4, 1, 2, 1, 1)

I am using the following command...

function [c25, c14, c2, vode, chi] = pde1dm_PS_OCP_v1(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_r)
% Piecewise PDE Implementation 
% Pulsing Simulation 
% Constants Used in the Nernst Equation   
n           = 2;                                             % No. of Electrons Involved
F           = 96485.3329;                                    % C/mol
R           = 8.314;                                         % J/(mol.K) or CV/(mol.K)
Temp        = 273.15+20;                                     % K
N           = 51; 
m           = 0; 
Area        = 0.0011;                                                     % Geometric Area (cm^2)
m_layer     = 4.13E-05;                                                   % m_layer (mol/cm^3)     
D_m         = 8.14E-08;                                                   % cm^2/s (De of polymer)
kk          = n*F*Area*l*m_layer; 
% DIMENSIONLESS PARAMETERS
chi         = linspace(0, 1, N);
tau         = linspace(0, tau_M1, N);                               % Dimensionless Time 
% FIRST POTENTIAL STEP    
E0          = 0.22;                                                 
E_set       = .350;                                      
epsilon1    = ((n*F)/(R*Temp))*(E_set-E0);                          % Dimensionless Potential 
c_M_ox      = 1/(1+exp(-epsilon1));                                 % Mox BC 
ic_arg_1    = {@(x)ones(size(N)) ; @(x)zeros(size(N))};
IC          =  @(x)PDE_ND_PS_EK_IC(x, ic_arg_1); 
BC          =  @(xl, yl, xr, yr, t)PDE_PS_EK_BC(xl, yl, xr, yr, c_M_ox); 
sol1        = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r), ... 
                    IC, BC, chi, tau);
 
c1          = sol1(:, :, 1);                                        % Substrate Conc.  
c2          = sol1(:, :, 2);                                        % Mox Conc.
               
% OPEN CIRCUIT POTENTIAL    
tau2         = linspace(tau_M1, tau_M2, N);                         % Dimensionless Time 
ic_arg_1     = {@(x)interp1(chi, sol1(N, :, 1), x) ; @(x)interp1(chi, sol1(N, :, 2), x)};
IC2          =  @(x)PDE_ND_PS_EK_IC(x, ic_arg_1); 
BC2          =  @(xl, yl, xr, yr, t, v)PDE_PSw_EK_BC_2(xl, yl, xr, yr, v); 
xOde         =  0;
ode_IC       =  @() ode_ICFun(tau_M1);
opts.vectorized='on';
% sol2         = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(x, t, c, DcDx, kappa, eta, gamma, mu, D_r), ... 
%                     IC2, BC2, chi, tau2);
[sol2, vode] = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(x, t, c, DcDx, kappa, eta, gamma, mu, D_r), ... 
               IC2, @PDE_PSw_EK_BC_2, chi, tau2, @odeFunc, ode_IC, chi(5));
               
% Concentration Profiles c(i, j, k)(Solutions) 
c14          = [c1; sol2(:, :, 1)];                                 % Substrate Conc.  
c25          = [c2; sol2(:, :, 2)];                                 % Mox Conc. 
c36          = 1-c25;                                               % Mred Conc. 
mox           = abs(sol2(:, 1, 2));
mred          = 1-mox;
E_PS          = E_set.*ones(1,N); 
for counter2  = 1:N
    E_OCP(counter2) = E0 + (((R*Temp)/(n*F).*(log(mox(counter2,:)./mred(counter2,:)))));
end
E_array       = [E_PS, E_OCP]; 
% This is the case finder function
% Video(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_m, N, chi, c14, c25, c36, E_array)
function [cc, ff, ss] = PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r)
 
    %                     S;          Mox;   
    cc              = [ D_r;           1]; 
    
    ff              = [ 1;             1].*DcDx;
    
    S_kin           = ((gamma/eta)*(kappa^2)*c(1)*c(2))./...
                      (gamma.*c(2).*(1+(mu.*c(1)))+c(1));
    M_kin           = ((kappa^2).*c(1).*c(2))./...
                      (gamma.*c(2).*(1+(mu.*c(1)))+c(1));
    ss              = [-S_kin;  -M_kin]; 
           
end
function c0     = PDE_ND_PS_EK_IC(x, ic_arg_1)
   
    %                         Substrate;                Mox;  
    c0              = [ic_arg_1{1}(x).';  ic_arg_1{2}(x).'];
end
% % % % % % % % % % % % % % % % For 0<tau_m<tau_1
function [pl, ql, pr, qr] = PDE_PS_EK_BC(xl, yl, xr, yr, c_M_ox)
    %                Substrate;         Mediator;
    pl     = [0               ;    yl(2)-c_M_ox];     
    
    ql     = [1               ;               0];
    
    pr     = [yr(1)-1         ;               0];
    
    qr     = [0               ;               1]; 
end
% % % % % % % % % % % % % % % % For tau_1<tau_m<tau_2
function [pl, ql, pr, qr] = PDE_PSw_EK_BC_2(xl, yl, xr, yr, t, v, vdot) 
    pl     = [0               ;               yl(2)-v];     
    
    ql     = [1               ;               0];
    
    pr     = [yr(1)-1         ;               0];
    
    qr     = [0               ;               1]; 
end
function f = odeFunc(t, v, vdot, x, c, DcDx, kk)
%     f = (DcDx*(1E03./(v(1)*(1-v(1)))))-vdot(1);
    f=DcDx(2) - 8/(c(2)*(1-c(2)))*vdot(1);
end
function v0=ode_ICFun(tau_0)
    v0 = ones(tau_0);
end
end

Hashim 2024-8-15，20:52

在 MATLAB Online 中打开

Hi Torsten,

I tried it with xOde=0 then I get this error. So I though might as well get a solution in the space just next to the surface at x=0 or chi=0. So I can run at chi(3)=0 or chi(2)=0 though it takes longer.

Warning: decic: Failed to obtain a converged set of consistent initial conditions. This might cause the ODE to DAE solver to fail in the first step. 
> In decicShampine (line 213)
In PDE1dImpl/solveTransient (line 212)
In pde1dm (line 123)
In pde1dm_PS_OCP_v1 (line 53) 
Warning: Failure at t=1.000000e+00.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-15) at time t. 
> In ode15i (line 392)
In PDE1dImpl/solveTransient (line 253)
In pde1dm (line 123)
In pde1dm_PS_OCP_v1 (line 53) 
Index in position 1 exceeds array bounds. Index must not exceed 1.
Error in pde1dm_PS_OCP_v1 (line 65)
    E_OCP(counter2) = E0 + (((R*Temp)/(n*F).*(log(mox(counter2,:)./mred(counter2,:)))));

Torsten 2024-8-15，22:12

编辑：Torsten 2024-8-15，22:12

The reason for your problems is obvious: c2(N,1) = 1.

That means that you divide by 0 in the expression f=DcDx(2) - 8/(c(2)*(1-c(2)))*vdot(1) .

In my opinion, there is something wrong in the mathematical formulation of your problem or in the choice of the starting values of c2 from phase 1 in phase 2.

请先登录，再进行评论。

Answer 2

Bill Greene 2024-8-20，13:30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2145369-how-to-input-a-novel-boundary-condition-for-a-coupled-pde-system#answer_1501794

编辑：Bill Greene 2024-8-20，13:53

在 MATLAB Online 中打开

Here is my code to solve this problem with pde1dm based on my understanding of the description and included code. I'm not at all sure I understand the problem definition and since I know nothing about the physics of the problem, I have no idea whether this solution is actually correct.

function matlabAnswers_8_15_2024
kappa=1; eta=1; gamma=1; mu=1; tau_M1=1; tau_M2=2; l=1; D_r=1;
pde1dm_PS_OCP_v1(kappa, eta, gamma, mu, tau_M1, tau_M2, 1, 1);
end
function [c25, c14, mox, vode] = pde1dm_PS_OCP_v1(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_r)
% Piecewise PDE Implementation 
% Pulsing Simulation 
% Constants Used in the Nernst Equation   
n           = 2;                                             % No. of Electrons Involved
F           = 96485.3329;                                    % C/mol
R           = 8.314;                                         % J/(mol.K) or CV/(mol.K)
Temp        = 273.15+20;                                     % K
Nx           = 21;
Nt=15;
m           = 0; 
Area        = 0.0011;                                                     % Geometric Area (cm^2)
m_layer     = 4.13E-05;                                                   % m_layer (mol/cm^3)     
D_m         = 8.14E-08;                                                   % cm^2/s (De of polymer)
k           = n*F*Area*l*m_layer;
% DIMENSIONLESS PARAMETERS
chi         = linspace(0, 1, Nx);
tau1         = linspace(0, tau_M1, Nt);                               % Dimensionless Time 
tau2         = [linspace(tau_M1+1e5*eps, tau_M2, Nt)];  
% FIRST POTENTIAL STEP    
E0          = 0.22;                                                 
E_set       = .350;                                      
epsilon1    = ((n*F)/(R*Temp))*(E_set-E0);                          % Dimensionless Potential 
c_M_ox      = 1/(1+exp(-epsilon1));                                 % Mox BC 
IC          =  @(x) PDE_ND_PS_EK_IC(x); 
BC          =  @(xl, yl, xr, yr, t,v,vDot) PDE_PS_EK_BC(xl, yl, xr, yr, t,v,vDot, c_M_ox, tau_M1, k); 
pdef= @(x, t, c, DcDx) PDE_ND_PS_EK(chi, tau1, c, DcDx, kappa, eta, gamma, mu, D_r);
xOde=0;
odeIcF= @() ode_ICFun(IC, chi);
[sol1,sol1Ode] = pde1dm(m, pdef, IC, BC, chi, tau1, @odeFunc, odeIcF, xOde);
%   figure; plot(tau1, sol1Ode, tau1, sol1(:,1,2));
%   title 'mox at left end as a function of time'; grid;
IC2=@(x) icFun2(x, chi, sol1);
odeIcF2= @() ode_ICFun(IC2, chi);
[sol2,sol2Ode] = pde1dm(m, pdef, IC2, BC, chi, tau2, @odeFunc, odeIcF2, xOde);
sol=[sol1; sol2];
solOde=[sol1Ode; sol2Ode];
tau=[tau1 tau2];
s          = sol(:, :, 1);                                        % Substrate Conc.  
mox          = sol(:, :, 2);                                        % Mox Conc.
figure; plot(tau, mox(:,1));
title 'mox at left end as a function of time'; grid;
figure; plot(tau, mox(:,end)); grid;
title 'mox at right end as a function of time';
figure; plot(chi, s(end,:));
title 's at final time'; grid;
figure; plot(chi, mox(end,:));
title 'mox at final time'; grid;
end
function [cc, ff, ss] = PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r)
%                     S;          Mox;
cc              = [ D_r;           1];
ff              = [ 1;             1].*DcDx;
S_kin           = ((gamma/eta)*(kappa^2)*c(1)*c(2))./...
    (gamma.*c(2).*(1+(mu.*c(1)))+c(1));
M_kin           = ((kappa^2).*c(1).*c(2))./...
    (gamma.*c(2).*(1+(mu.*c(1)))+c(1));
ss              = [-S_kin;  -M_kin];
end
function c0     = PDE_ND_PS_EK_IC(x)
%                         Substrate;                Mox;  
c0              = [1 0]';
end
function c0 = icFun2(x, chi, sol1)
s1Final=squeeze(sol1(end,:,:));
c0=interp1(chi, s1Final, x)';
end
function [pl, ql, pr, qr] = PDE_PS_EK_BC(xl, yl, xr, yr, t,v,vDot, c_M_ox, tau_M1, k)
% ---------------------
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% |      |            |
% -------|-------------
%        pl           pr
%                Substrate;         Mediator;
if t <= tau_M1
    pl     = [0               ;    yl(2)-c_M_ox];
    ql     = [1               ;               0];
else
    moxR=yr(2);
    dMoxDTau=vDot;
    p=-k/(moxR*(1-moxR))*dMoxDTau;
    pl     = [0               ;    p];
    ql     = [1               ;    1];
end
pr     = [yr(1)-1         ;               0];
qr     = [0               ;               1];
end
function F=odeFunc(t,v,vdot,x,u)
F=v-u(2);
end
function v0=ode_ICFun(icf,chi)
x0=chi(1);
c0=icf(x0);
v0 = c0(2);
end