Boundary condtions for an index reduced DAE system using ode solvers

Hey,

so i have the following set of differential equations. Its called the poisson nerst planck equation and i am trying to solve it using the ode solvers of matlab.

with the boundary conditions for are neuman boudnary conditions

(it depends on i)

and for ϕ are dirichlet boundary conditions

Now I have allready tried to formulate this problem into a system of frist oder differential equation as follows:

Now i insert and in and solve for

The complete ode system should consist of these 3 first oder equations now:

with boundary conditions:

as you can see I have only 3 first oder equations but 4 boundary conditions. aditionally i dont have boundary conditions for for each but rather two boundary conditions for and two boudnary condition for .

Now i am clueless on how to implement the boundary conditions for this system of equations with a suitible ode solver for matlab. Is there a way of implementing two dirichlet and two neuman boundary conditions?

As an extra. I want to solve this set of equations for 3 adjacient regions within the overall domain of as , where for and and for . This gives the problem a stiff character becasue almost instantanious increas of at these internal regions.

[Torsten]

Why did you proceed differently than suggested to your previous question ?

[Berat Cagan Türkmen]

Hey,

so i tried you suggestion but i can not find physically corrent boundary conditions. As i have 2 first oder equations in 3 regions i need 6 boundary conditions. I was able to set the boundary conditions at the outermost boundaries. And i also set the countinuity boudnary for the conectration c_i or y(1). And this also gave me the jumps i wanted to see.You can also see in the screenshot provieded.

but i still need 2 boudnary conditions for dc_i/dx or y(2) which i dont know how to come up with. I know that the gradient of c_i has change at the interface but as you can see form the screenshot its just stays the same.

now i am trying a different approach were i dont "imposed" discontinuity. But this gives me this versty stiff problem which i initially wanted to avoid.

if oyu are interested you can read more about it here:

https://www.comsol.com/blogs/how-to-model-ion-exchange-membranes-and-donnan-potentials/

%parameters
DH = 3.33*10^-9;
DHM = 3.5*10^-10;
F = 96485;
I = 100;
R = 8.314;
T = 300;
f = F/(R*T);
zH = 1;
RH =I/(zH*F*DH);
PHI = 10;
PHIM = 1;
p = [DH DHM F I R T f zH RH PHI PHIM];
%mesh
N = 100;    %step size
L= 0.005;   %chanel width in m
dx = L/N;   %step size in m
xmesh1 = [0:dx:L];
xmesh2 = [L:dx:2*L];
xmesh3 = [2*L:dx:3*L];
xmesh = [xmesh1,xmesh2,xmesh3];
%initial condition
yinit = [10; 10];
sol = bvpinit(xmesh,yinit);
%solver
sol = bvp5c(@(x,y,r) fun(x,y,r,p), @bc, sol);
%plot
plot(sol.x,sol.y(1,:))
ylabel('y1(x)')

function dydx = fun(x,y,region,p) % equations being solved
%n = p(1);
%eta = p(4);
DH = p(1);
DHM = p(2);
F = p(3);
I = p(4);
R = p(5);
T = p(6);
f = p(7);
zH = p(8);
RH = p(9);
PHI = p(10);
PHIM = p(11);

dydx = zeros(2,1);

switch region
    case 1    % x in [0 1]
        dydx(1)=y(2);
        dydx(2) = -zH*f*y(2)*PHI;
    case 2    % x in [1 lambda]
        dydx(1)=y(2);
        dydx(2) = -zH*f*y(2)*PHIM;
    case 3    % x in [1 lambda]
        dydx(1)=y(2);
        dydx(2) = -zH*f*y(2)*PHI;
end
end
%-------------------------------------------
function res = bc(YL,YR) % boundary conditions
DH = 3.33*10^-9;
DHM = 3.5*10^-10;
F = 96485;
I = 100;
R = 8.314;
T = 300;
f = F/(R*T);
zH = 1;
RH =I/(zH*F*DH);
PHI = 10;
PHIM = 1;
% res = [YL(2,1)+RH   
%        YR(1,1) - YL(1,2)+100%DH*YR(2,1)+zH*DH*f*YR(1,1)*PHI - DHM*YL(2,2)-zH*DHM*f*YL(1,2)*PHIM
%        YR(2,1) - YL(2,2)%DH*YR(2,1)+zH*DH*f*YR(1,1)*PHI - DHM*YL(2,2)-zH*DHM*f*YL(1,2)*PHIM%
%        YR(1,2) - YL(1,3)-100%DHM*YR(2,2)+zH*DHM*f*YR(1,2)*PHIM - DH*YL(2,3)-zH*DH*f*YL(1,3)*PHI
%        YR(2,2) - YL(2,3)%DHM*YR(2,2)+zH*DHM*f*YR(1,2)*PHIM - DH*YL(2,3)-zH*DH*f*YL(1,3)*PHI% 
%        YR(2,3)]
res = [YL(2,1)+RH   
       (DH*YR(2,1)+zH*DH*f*YR(1,1)*PHI) + (DHM*YL(2,2)+zH*DHM*f*YL(1,2)*PHIM)
       YR(2,1) - YL(2,2)%(DH*YR(2,1)+zH*DH*f*YR(1,1)*PHI) - (DHM*YL(2,2)+zH*DHM*f*YL(1,2)*PHIM) %%HERE%%
       (DHM*YR(2,2)+zH*DHM*f*YR(1,2)*PHIM) + (DH*YL(2,3)+zH*DH*f*YL(1,3)*PHI)
       YR(2,2) - YL(2,3)%(DHM*YR(2,2)+zH*DHM*f*YR(1,2)*PHIM) - (DH*YL(2,3)+zH*DH*f*YL(1,3)*PHI)  %%HERE%%
       DH*YL(2,3)+zH*DH*f*YL(1,3)*PHI];      
end
%-------------------------------------------

[Torsten]

Both YR(1,3) and YR(2,3) are missing in your specification of boundary conditions in your above code. That is not possible - you have to set a condition at the outer right boundary.

[Berat Cagan Türkmen]

Thank you for the hint, i changes the code acording to your suggestion. But im still left with my initial problem, which is that i dont know the boundary of YR(2,1)/YL(2,2) and YR(2,2)/YL(2,3).

Additionally I am worried that this approach is the wrong way to handle this problem. because the jumping conentration at the boundary should be because of the which apprears in the membrane (center) region. but in this model i did not even implement that yet. Which is why i am almost certain that its the wrong approach.

That is why im trying a new appreaoch. I appreaciate any comments type of comments or advices.

%parameters
DH = 3.33*10^-9;
DHM = 3.5*10^-10;
F = 96485;
I = 100;
R = 8.314;
T = 300;
f = F/(R*T);
zH = 1;
RH =I/(zH*F*DH);
PHI = 10;
PHIM = 1;
p = [DH DHM F I R T f zH RH PHI PHIM];
%mesh
N = 100;    %step size
L= 0.005;   %chanel width in m
dx = L/N;   %step size in m
xmesh1 = [0:dx:L];
xmesh2 = [L:dx:2*L];
xmesh3 = [2*L:dx:3*L];
xmesh = [xmesh1,xmesh2,xmesh3];
%initial condition
yinit = [10; 10];
sol = bvpinit(xmesh,yinit);
%solver
sol = bvp5c(@(x,y,r) fun(x,y,r,p), @bc, sol);
%plot
plot(sol.x,sol.y(1,:))
ylabel('y1(x)')

function dydx = fun(x,y,region,p) % equations being solved
%n = p(1);
%eta = p(4);
DH = p(1);
DHM = p(2);
F = p(3);
I = p(4);
R = p(5);
T = p(6);
f = p(7);
zH = p(8);
RH = p(9);
PHI = p(10);
PHIM = p(11);

dydx = zeros(2,1);

switch region
    case 1    % x in [0 1]
        dydx(1)=y(2);
        dydx(2) = -zH*f*y(2)*PHI;
    case 2    % x in [1 lambda]
        dydx(1)=y(2);
        dydx(2) = -zH*f*y(2)*PHIM;
    case 3    % x in [1 lambda]
        dydx(1)=y(2);
        dydx(2) = -zH*f*y(2)*PHI;
end
end
%-------------------------------------------
function res = bc(YL,YR) % boundary conditions
DH = 3.33*10^-9;
DHM = 3.5*10^-10;
F = 96485;
I = 100;
R = 8.314;
T = 300;
f = F/(R*T);
zH = 1;
RH =I/(zH*F*DH);
PHI = 10;
PHIM = 1;

res = [YL(2,1)+RH   
       (DH*YR(2,1)+zH*DH*f*YR(1,1)*PHI) + (DHM*YL(2,2)+zH*DHM*f*YL(1,2)*PHIM)
       YR(2,1) - YL(2,2)%HERE
       (DHM*YR(2,2)+zH*DHM*f*YR(1,2)*PHIM) + (DH*YL(2,3)+zH*DH*f*YL(1,3)*PHI)
       YR(2,2) - YL(2,3)%HERE
       DH*YR(2,3)+zH*DH*f*YR(1,3)*PHI];      
end
%-------------------------------------------

[Torsten]

Well, you know best the physical background of the problem. I cannot help you in this respect.

But I cannot imagine that solving the additional equation for phi can help to circumvent specifying adequate transmission conditions for the concentrations.

By the way: In your problem formulation from above, you forgot the equation

dy4/dx = -F/eps*sum_(z_i*c_i) + rho_fix

So what is the problem ? Four first-order equations with four boundary conditions.

[Berat Cagan Türkmen]

to maybe further clearafy my question.

I am looking to implmeenting two neuman boundary condition for c_i. which means i need to define y(2) at 0 and at L.

But form what i understand the solver want me to give 1 boundary ocndition for y(1) and one condition for y(2).

[Torsten]

But form what i understand the solver want me to give 1 boundary ocndition for y(1) and one condition for y(2).

No. The number of boundary conditions must equal the number of equations. That's the only restriction for that the solver accepts your settings. Whether these settings make sense and give a unique solution is a different question.

Same for the transmission conditions.

But giving both conditions for the derivative is dangerous in general because (at least for linear ODEs) together with c, also c+constant is a solution.

[Berat Cagan Türkmen]

Thanks for the corrention on the ode system.

No. The number of boundary conditions must equal the number of equations. That's the only restriction for that the solver accepts your settings. Whether these settings make sense and give a unique solution is a different question.

Then i might not know the correct way of specifying boundary condition.

If you look at the code for the system you can see y0, which are supposed to be my boundary conditions.

here is how I understand it based on the documentation.

y0(1) is the boundary condition for y(1) which is c_i. Specifing a number here would therefore give me a dirichlet boundary condition.

for y0(2) which is dc/dx i would define a neuman boundary condition for c_i.

And i realize that i also dont a know how to specify a coundition at a certain boundary.(for 0 or for L)

x = [0; L]; % several periods

y0 = [0;0;0;0];

[t,y] = ode15s(@f,x,y0);

figure;

plot(t,y(:,1));

function dydt = f(t,y)

%parameters

DH = 3.33*10^-9;

DHM = 3.5*10^-10;

F = 96485;

I = 100;

R = 8.314;

T = 300;

f = F/(R*T);

zH = 1;

RH =I/(zH*F*DH);

PHI = 10;

PHIM = 1;

dydt=zeros(4,1)

if x<=L/3

zfix = 0;

cfix = 0;

dydt(1)=y(2);

dydt(2)=-zH*f*y(1)*y(4)+zH*f*(F/eps)*(zH*y(1))+zfix*cfix;

dydt(3)=y4;

dydt(4)=zH*f*(F/eps)*(zH*y(1))+zfix*cfix;

elseif (x >= L/3) && (x <= 2/3*L)

zfix = 2;

cfix = 2;

dydt(1)=y(2);

dydt(2)=-zH*f*y(1)*y(4)+zH*f*(F/eps)*(zH*y(1))+zfix*cfix;

dydt(3)=y4;

dydt(4)=zH*f*(F/eps)*(zH*y(1))+zfix*cfix;

elseif (x >= 2/3*L) && (x <= L)

zfix = 0;

cfix = 0;

dydt(1)=y(2);

dydt(2)=-zH*f*y(1)*y(4)+zH*f*(F/eps)*(zH*y(1))+zfix*cfix;

dydt(3)=y4;

dydt(4)=zH*f*(F/eps)*(zH*y(1))+zfix*cfix;

end

% -----------------------------------------------------------------------

[Torsten]

ode15s solves initial value problems, not boundary value problems.

You have a boundary value problem.

That's why you have to use bvp4c or bvp5c.

Or program your own code.

[Berat Cagan Türkmen]

%parameters

DH = 3.33*10^-9;

DHM = 3.5*10^-10;

F = 96485;

I = 100;

R = 8.314;

T = 300;

f = F/(R*T);

zH = 1;

RH =I/(zH*F*DH);

PHI = 10;

PHIM = 1;

p = [DH DHM F I R T f zH RH PHI PHIM];

%mesh

N = 100; %step size

L= 0.005; %chanel width in m

dx = L/N; %step size in m

xmesh = [0:dx:3*L];

%initial condition

yinit = [100;1;0;0];

sol = bvpinit(xmesh,yinit);

%solver

%sol = bvp5c(@(x,y,r) fun(x,y), @bc, sol);

eps = 1;

for i = 2:13

eps = eps/10;

sol = bvp5c(@(x,y,r) fun(x,y,eps), @bc, sol);

end

%plot

plot(sol.x,sol.y(1,:),sol.x,sol.y(3,:));

ylabel('y1(x)')

function dydt = fun(x,y,eps)

%parameters

DH = 3.33*10^-9;

DHM = 3.5*10^-10;

F = 96485;

I = 0.1;

R = 8.314;

T = 300;

f = F/(R*T);

zH = 1;

RH =I/(zH*F*DH);

PHI = 10;

PHIM = 1;

%eps = 1e1;

L=0.015;

dydt=zeros(4,1);

if x<=L/3

zfix = 0;

cfix = 0;

dydt(1)=y(2);

dydt(2)=-zH*f*y(1)*y(4)+zH*f*(F/eps)*(zH*y(1))+zfix*cfix;

dydt(3)=y(4);

dydt(4)=zH*f*(F/eps)*(zH*y(1))+zfix*cfix;

elseif (x >= L/3) && (x <= 2/3*L)

zfix = 1;

cfix = 2;

dydt(1)=y(2);

dydt(2)=-zH*f*y(1)*y(4)+zH*f*(F/eps)*(zH*y(1))+zfix*cfix;

dydt(3)=y(4);

dydt(4)=zH*f*(F/eps)*(zH*y(1))+zfix*cfix;

elseif (x >= 2/3*L) && (x <= L)

zfix = 0;

cfix = 0;

dydt(1)=y(2);

dydt(2)=-zH*f*y(1)*y(4)+zH*f*(F/eps)*(zH*y(1))+zfix*cfix;

dydt(3)=y(4);

dydt(4)=zH*f*(F/eps)*(zH*y(1))+zfix*cfix;

end

end

function res = bc(YL,YR) % boundary conditions

DH = 3.33*10^-9;

DHM = 3.5*10^-10;

F = 96485;

I = 10;

R = 8.314;

T = 300;

f = F/(R*T);

zH = 1;

RH =I/(zH*F*DH);

PHI = 10;

PHIM = 1;

res = [YL(2)+RH

YR(2)

YL(3)

YR(3)-4

];

end

%-------------------------------------------

so i tried to implemnt the problmem but i am not getting satisfying results.

the variable eps is around 10^-12 which is suposed to cause an almost instantanious increase at the internal interfaces. (like the jump condition from my other post). Due to the high eps value the problem becomes stiff.

I tried to iteratively increase the value of eps as shown in this docuentation:

https://de.mathworks.com/help/matlab/math/solve-bvp-using-continuation.html

but with no success.

Are there any other ways i can improve the solvability?