How can I make a bvp4c solution more stable?

Hello community,

I have a problem with the stability of the BVP4C-solution when the derivative becomes very small. The equation system describes the molar flow of four components along a membrane. The gradient of each species along the membrane is the product of the respective permeability Q and the difference of partial pressure between the high pressure side (retentate) and low pressure side (permeate). Problems occur when the partial pressure of the fast diffusing component on both sides is more or less equal.

I am solving a boundary value problem of the following form

y'= dy/dx = f(y,yo)

y is a vector of the molar flow of 4 species in the retentate side of the membrane. At x=0 I have to solve an implicit equation system for a set of parameters yPo

g(yo,yPo)=0

and calculate the derivative

y'(x=0)=f(yo,yPo)

yo is y(x=0), an unknown parameter (my left boundary value - the retentate molar flows at the end of the membrane) which I fit with the bvp4c-solver. I know the right boundary value at x=1 where the feed enters the membrane on the retentate side.

The ODE is quite simple:

function dydx=ydot(x,y,yo,n)
  Q = [2.3369E-08; 1.17095E-06; 3.81087E-08; 2.36619E-07]*n;
  pR = 1e6;
  pP = 1e5;
  if x==0
      options=optimset('display','none');
      yPo=fsolve(@(yPo) PermOmega(yPo),Q/sum(Q),options);
      dydx=Q.*(y/sum(y)*pR-yPo*pP);
  else
      dydx=Q.*(y/sum(y)*pR-(y-yo)/sum(y-yo)*pP);
  end

    function g=PermOmega(yPo)
        %Permeate composition at x = 0
        g=yPo-Q.*(yo/sum(yo)*pR-yPo*pP)/sum(Q.*(yo/sum(yo)*pR-yPo*pP));
    end
end

My initial guess for the solution is a straight line between the two boundary values with a mesh-size of 10. The problem is in the second line of the code, where the permeability is multiplied with the number of membranes n. The solution is stable for n=1 (one membrane) but when I increase the number of membranes the partial pressure become similar and the ODE becomes unstable (see pictures)

What can I do to improve stability when the gradient over the membrane becomes small?

So far I have provided the Jacobian matrices (using implicit differentiation at x=0), which made it faster but not more accurate. I have changed the relative and absolute error tolerances of the BVP4C-solver. I have tried BVP5C. And I have boiled down my original ODE-system with three stages each with 7 variables (pressure in retentate and permeate plus retentate temperature) to the essential system with only 4 variables.

The reason why I need stability is that I want to use the solution in a FMINCON-optimizer which will sometimes select parameters which will lead to unstable solutions.

Any ideas?

Answer 1

[Torsten]

Thanks for the detailed explanation of the problem.

I used the following problem formulation in COMSOL Multiphysics and it works without problems:

Solution variables:

y, yo

Equations:

dydx=Q.*(y/sum(y)*pR-(y-yo)/sum(y-yo)*pP)
dyodx = 0

Initial conditions:

y=yo0+x*(yalpha-yo0) 
yo=yo0

Boundary conditions:

yo(x=0)=y
y(x=1)=yalpha

Note that I solve for 8 unknown functions, not only 4.

Note further that ydot is only called at the inner grid points, not at x=0. Thus a division by zero does not occur.

Best wishes

Torsten.

Answer 2

[Torsten]

You must specify the retentate molar flows at x=0 in the boundary routine of bvp4c, not in the routine where you specify the ODE.

Best wishes

Torsten.

Answer 3

[JK]

Hi Torsten,

thanks for the answer. Here is the code for calling the BVP4C:

% initial conditions at x=1 (retentate inlet)
yalpha = [0.3;0.25;0.01;0.01]; 

% estimated retentate values at the end of the membrane (x=0)
yo0=[0.95;0.60;0.92;0.70].*yalpha;

% Boundary value solver
FJac=@(x,y,yo) Jacfun(x,y,yo,n);  
BCJac={[eye(4);zeros(4)],[zeros(4);eye(4)],[-eye(4);zeros(4)]};
bcfun = @(ya,yb,yo) [ya-yo; yb-yalpha];
initguess = @(x) yo0+x*(yalpha-yo0);
odefun =@(x,y,yo) ydot(x,y,yo,n);
solinit = bvpinit(linspace(0,1,11),initguess,yo0);
options = bvpset('stats','on','RelTol',1e-3,'FJacobian',FJac,'BCJacobian',BCJac); 
sol=bvp4c(odefun,bcfun,solinit,options);

Jacfun computes two Jacobian matrices dy'/dy and dy'/dyo (code not shown) each matrix has to be calculated separately for x=0 and x>0.

bcfun = @(ya,yb,yo) [ya-yo; yb-yalpha] computes the boundary conditions. yo is the parameter which is computed by the boundary value solver.

I specify y(x=0) in the boundary condition and I specify dy/dx(x=0) in the ODE-equation. I still do not know how to solve the issue. I think it has something to do with the stiffness of the ODE.

kind regards

JK

[Torsten]

Try

% initial conditions at x=1 (retentate inlet)
yalpha = [0.3;0.25;0.01;0.01]; 

% estimated retentate values at the end of the membrane (x=0)
yo0=[0.95;0.60;0.92;0.70].*yalpha;

% Boundary value solver
bcfun = @(ya,yb,yo) bc(ya,yb,yo,n,yalpha); 
initguess = @(x) yo0+x*(yalpha-yo0);
odefun =@(x,y,yo) ydot(x,y,yo,n);
solinit = bvpinit(linspace(0,1,11),initguess,yo0);
options = bvpset('stats','on','RelTol',1e-3); 
sol=bvp4c(odefun,bcfun,solinit,options);

function dydx=ydot(x,y,yo,n)
  Q = [2.3369E-08; 1.17095E-06; 3.81087E-08; 2.36619E-07]*n;
  pR = 1e6;
  pP = 1e5;
  dydx=Q.*(y/sum(y)*pR-(y-yo)/sum(y-yo)*pP);
end

function res=bc(ya,yb,yo,n,yalpha)
  Q = [2.3369E-08; 1.17095E-06; 3.81087E-08; 2.36619E-07]*n;
  pR = 1e6;
  pP = 1e5;
  options=optimset('display','none');
  yPo=fsolve(@(yPo) PermOmega(yPo),Q/sum(Q),options);
  res = [ya-yPo;yb-yalpha];

      function g=PermOmega(yPo)
          %Permeate composition at x = 0
          g=yPo-Q.*(yo/sum(yo)*pR-yPo*pP)/sum(Q.*(yo/sum(yo)*pR-yPo*pP));
      end
  end

Answer 4

[JK]

Hi Torsten,

thank you for the time you spend on this problem. The whole issue is about the calculation of yo, the retentate flow at x=0 with the boundary value solver. yPo is the composition of the permeate at x=0, its unit is mol/mol, whereas the flow in the retentate y is in mol/s. I thought that might be confusing. Historically y is used in membrane science for the permeate composition, in thermodynamics for gas composition in general and in math for the dependent variable.

The ODE descibes the flow over the membrane, it is permability (Q) times partial pressure difference over the membrane (pressure retentae times molar fraction in retentate minus pressure permeate times molar fraction in permeate). the molar fraction in the retentate yR is

yR= nR/sum(nR) (nR is called y in my program - that is confusing, I admit)

the molar flow in the permeate at x can be calulated by balancing from x=0 to x

nP = nR-nRo

so

yP= (nR-nRo)/sum(nr-nRo) (or (y-yRo)/sum(y-yRo) in the code)

at x=0 nr=nro, the permeate flow is zero and the transportation eqation becomes 0/0 NaN. So at x=0 I need a special equation to avoid the NaN. At x=0 the permeate flow is zero, but the permeate still has a composition. It is

yPo = n*/sum(n*)

n* is the flow over the membrane and it is

n* = Q*(pR yo/sum(yo) - pP yPo)

and this is calculated with the solver routine PermOmega.

If I could calculate the retentate concentration at x=0 (yo) directly, I would not have to solve the BV-problem. There might be an algebraic solution for this but when I add the pressure and temperautre profiles to the mix, I will need to solve the ODE.

So how can I improve the stability of the bvp4c-routine? Yesterday I tried to gradually increase the number of membranes and used the old solution as a new starting point (Continuation). It didn't work. When I increase the number of membranes from 1.55 to 1.56 the system drifts into chaos and the smooth profile will becomme a chaotic one. I also changed the mesh size. A smaller mesh size with more points at x close to 1 will actually improve stability a bit.

Any more ideas how to make this work? It must be possible somehow.

Answer 5

[JK]

ok, can't wait to try it this evening, will keep you updated.

EDIT: I get a singular jacobian when the solver reaches x=1. In my boundary-conditions I do not define the right boundary value of yo. Interestingly I can only pass 8 variables into the Boundary functions.

I define the new y as [y_old;yo], so the new y = y(1:4) and yo = y(5:8).

here is the code, it became quite simple:

clear
global yalpha
yalpha = [0.3;0.25;0.01;0.01]; % inlet molar flow
n=1; % no of membranes
Q = [2.3369E-08; 1.17095E-06; 3.81087E-08; 2.36619E-07]*n; % Permeability
yo0=[0.95;0.60;0.92;0.70].*yalpha; % estimated values for outlet flow

% Boundary value solver

bcfun = @(ya,yb) bc(ya,yb);
initguess = @(x) init(x,yo0);
odefun =@(x,y,yo) ydot(x,y,n);
solinit = bvpinit(linspace(0,1,11),initguess);
options = bvpset('stats','on','RelTol',1e-4);

sol=bvp4c(odefun,bcfun,solinit,options);

% functions -------------
function IG=init(x,yo0)   % initial guess of mesh
    global yalpha
    IG=zeros(8,1);
    IG(1:4)=yo0+x*(yalpha-yo0);
    IG(5:8)=yo0;
end

function res=bc(ya,yb)     % boundary value function
   global yalpha
   res=zeros(8,1);
   res(1:4)=ya(1:4)-ya(5:8);    %y at left BV (x=0) y=yo
   res(5:8)=yb(1:4)-yalpha;     %y at right BV (x=1) y=yalpha
end

function dydx=ydot(x,y,n)       % derivative
    dydx=zeros(8,1);
    Q = [2.3369E-08; 1.17095E-06; 3.81087E-08; 2.36619E-07]*n;
    pR = 1e6;
    pP = 1e5;
    dydx(1:4)=Q.*(y(1:4)/sum(y(1:4))*pR-(y(1:4)-y(5:8))/sum(y(1:4)-y(5:8))*pP);
    dydx(5:8)=zeros(4,1);
end

I think we're almost there...

[Torsten]

Use

odefun =@(x,y) ydot(x,y,n);

instead of

odefun =@(x,y,yo) ydot(x,y,n);

Best wishes

Torsten.

[JK]

still singular Jacobian. yb(5:8) is not defined in boundary-function. The dres/dyb(5:8) lines in the jacobian of the boundary-values are empty.

Answer 6

[JK]

tada:

The problem was the solver for yPo at x=0. when n>>1 sometimes it comes up with solutions where min(yPo)<0 or max(yPo)>1 which makes no sense phyiscally. So I calcualte yPo for n = 1 or whatever is a stable solution and then use it as a new startingpoint for future calculations of yPo. I use continuation to extend the solution to the desired number of membranes.

Here is the code:

yalpha = [0.3;0.25;0.01;0.01]; % input molar flow 

N=5; % number of membranes

global yPo0
Q = [2.3369E-08; 1.17095E-06; 3.81087E-08; 2.36619E-07]; permeabilities
yPo0=Q/sum(Q); initial guess

% estimated values at the end of the membrane (x=0)
yo0=[0.95;0.60;0.92;0.70].*yalpha;

% start
% ---------------------------------------------------------------------
tic;

% Boundary value solver

xx=linspace(0,1,10);
xx=xx.^0.5;         % moves mesh points towards x =1

figure
BCJac={[diag(ones(1,4));zeros(4,4)],[zeros(4,4);diag(ones(1,4))],[diag(-ones(1,4));zeros(4,4)]};
bcfun = @(ya,yb,yo) [ya-yo; yb-yalpha];
initguess = @(x) IG(x,yo0,yalpha);
sol = bvpinit(xx,initguess,yo0);
nn=(N-1)/10;
for n = 1:nn:N
    % Permeability
    Q = [2.3369E-08; 1.17095E-06; 3.81087E-08; 2.36619E-07]*n;
    FJac=@(x,y,yo) Jacfun(x,y,yo,n);
    odefun =@(x,y,yo) ydot(x,y,yo,n);
    options = bvpset('stats','on','RelTol',1e-3,'BCJacobian',BCJac); %'FJacobian',FJac,
    sol=bvp4c(odefun,bcfun,sol,options);
    yy=deval(sol,xx);
    plot(xx,yy(2,:));
    axis([0,1,0,0.3]);
    title(num2str(n));
    pause(0.2);
    sol=bvpinit(xx,@(x) deval(sol,x),sol.parameters); % makes sure that the mesh is not bloating with mesh points as the calculation goes on
end

disp(sol);
toc;
x=sol.x;
y=sol.y;

plot(x,y);
title([num2str(n),' membranes']);
xlabel('dimensionless length');
ylabel('mol flow retentate');

function dydx=ydot(x,y,yo,n)
    Q = [2.3369E-08; 1.17095E-06; 3.81087E-08; 2.36619E-07]*n;
    pR = 1e6;
    pP = 1e5;
    if x==0
        global yPo0
        options=optimset('display','none');
        if min(yPo0)<0 || max(yPo0)>1
            yPo0=[0.2;0.75;0.03;0.02];
            beep; % if you hear this: not good
        end
        yPo=fsolve(@(yPo) PermOmega(yPo),yPo0,options);
        yPo0=yPo;   % that was the breakthrough !
        dydx=Q.*(y/sum(y)*pR-yPo*pP);
    else
        dydx=Q.*(y/sum(y)*pR-(y-yo)/sum(y-yo)*pP);
    end

      function f=PermOmega(yPo)
          %Permeate composition at x = 0
          f=yPo-Q.*(yo/sum(yo)*pR-yPo*pP)/sum(Q.*(yo/sum(yo)*pR-yPo*pP));
      end

end

   function [FJ,FJo]=Jacfun(x,y,yo,n)
      FJ=Jacy(x,y,yo,n);
      FJo=Jacyo(x,y,yo,n);  % external functions for Jacobians
   end

   function f=IG(x,yo0,yalpha)  % initial guess only used for the first calculation
      f= yo0+x*(yalpha-yo0);
      f(2)=yo0(2)+x^0.5*(yalpha(2)-yo0(2));  % takes the convex nature of y(2) into account
   end

It is funny that the first thing they tell you is that when numeric solutions are calculated it is all about the starting values. And you think you understand but after spending dozen of nights staring at computer screens you finally realize: it is all about starting values!

Thanks Torsten, You have been a great help!

Answer 7

[Nadeem Abbas]

Hello Dear

How to find the stability of the solutions of BVP4C method?

[Chandan K G]

by taking different initial guesses according to literature