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

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

Answer 1

[Torsten]

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

[Hashim]

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

[Torsten]

Your comment doesn't change my answer given.

I understood everything as you describe it.

The defined ode in my answer solves for m_ox in chi = 0. I guess dm_ox/dtau in your boundary condition means dm_ox/dtau in chi = 0 and m_ox(1) means m_ox at chi = 1.

Of course, you cannot set m_ox = 1 at chi = 1 and t = 0 as initial condition because this will immediately give a division by zero in the boundary condition.

[Bill Greene]

I took a quick look at your equations and code.

Do I understand correctly that the only thing that changes in the

two time spans is the LHS BC for m_0x? If so why do you call pde1dm twice?

Your approach will probably work but seems unnecessarily complex to me.

To implement your BC, I believe you will need ODE coupling points at

zero and one. I have no idea why you have one at .9.

[Torsten]

@Bill Greene

I couldn't find the ordering in your documentation. Can you clarify ?

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

Does your pde1dm work with full or banded Jacobian ? I ask this because the new boundary condition couples values of c in x = 0 and x = 1 which would require a full Jacobian.

[Hashim]

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]

Hi @Bill Greene yes it is the boundary condition that changes but I also need to apply the initial condition using the concentration profile from the last time step. What would be the proper way to go about this then in your opinion any examples I can follow?

[Bill Greene]

Yes, it looks like that is one of the many areas where the pde1dm doc could

use some improvements!

But, in the meantime, your assumption is correct:

% Assuming first index in DcDx and c is equation, second index is coupling point in xOde vector

pde1dm provides a sparsity pattern to ode15s that assumes all ODE equations are fully

coupled to all of the discretized PDE equations. But it does not try to analyze the ODE

equations to recognize that they cause other terms in the PDE equations to be coupled; so,

in that sense, the jacobian is only approximate. There is no such approximation in the

calculation of residual of the system so, if at a time step, the newton iteration converges,

it is correct.

I have noticed in these implicit ODE solvers that it is often possible for the solution to

converge with a fairly rough, approximate jacobian. But I strongly suspect it is possible to

construct test cases that fail to converge. And in most cases, the approximate jacobian

likely slows down the rate of convergence.

[Torsten]

@Hashim

Why don't you define the ODE as below and not as I wrote ?

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

[Hashim]

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

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>

<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>

<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>

<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>

<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

[Torsten]

It seems that you cannot set an ordinary differential equation that uses values from your PDE solution in two different points simultaneously. "odeFunc" is called for each coupling point separately (see above). Or there is a setup error elsewhere in your code.

But maybe m_ox at chi = 1 remains constant during the integration and you can insert it directly for m_ox(1) in the ODE equation. Then you only have one coupling point, namely chi = 0 (xOde = 0), and you can proceed as suggested.

I just see that you changed your boundary condition in your first message. Why don't you inform us about that ? Now the setup is quite simple as you only have chi = 0 as coupling point.

[Hashim]

Hi, my bad I thought i'd said it in my earlier replies but I guess it did not get updated or something. The mistake was there because i have another script with the method implemented in MOL. That script also needs just the BC fixed.

Anyways I tried your odeFunc and set xOde=0 but it did not work for me. Did you change anything else as well? Or perhaps it is just not possible.

[Bill Greene]

Your ODE definition function is incorrectly defined.

Look at section 3.2.1 in the manual.

[Hashim]

My God good Sirs, it seems to be working. Here is the updated odeFunc function. I will play with the script a bit more so hopefully you will bother with a few more q's.

I have a question though why can't I define the vdot as cdot. I mean on one hand I am using c's as my original PDE but then I have to define a separate variable v for my ODE and solve for it. It's a bit weird in my opinion. And why is vode a matrix instead of an array.

function f = odeFunc(t, v, vdot, x, c, DcDx)

%     f = (DcDx*(1E03./(v(1)*(1-v(1)))))-vdot(1);
x
size(DcDx)
size(c)
size(vdot)
    f=DcDx(2) - 10/(c(2)*(1-c(2)))*vdot(1);

end

[Torsten]

I changed your code from above and ran it. I'm surprised that v remains at its initial value.

I have a question though why can't I define the vdot as cdot. I mean on one hand I am using c's as my original PDE but then I have to define a separate variable v for my ODE and solve for it. It's a bit weird in my opinion.

The possible boundary conditions for c are of the form p + q*f = 0. There is no way to define a differential equation in time for c as a Dirichlet boundary condition. So we must artificially introduce a variable v (which is c at chi = 0) and solve an ODE for it in order to set a boundary condition for c of the above form (p = c - v, q = 0).

And why is vode a matrix instead of an array.

I get an 21x1 array for v (see above).

[Hashim]

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

[Torsten]

At least these two lines are wrong in your code:

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

instead of

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

and

[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);

instead of

[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));

[Hashim]

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]

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]

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

[Torsten]

OP changed his problem description from

moxR=yr(2);

to

moxL=yl(2);