How can I introduce constraints to the nature of solutions provided by a boundary value problem solver such as bvp4c or bvp5c?

I have two coupled differential equations that I solve using bvp4c. The equations relate to concentrations of two different constituents C_Ga and C_As:

The values of all other equation parameters (D_Ga, D_As etc) are known. For certain values of the equation parameters, the solution produced by bvp4c is perfect. However, for other parameter values (that are nonetheless physically reasonable), the solution produced by bvp4c is clearly wrong. For example, the code produces negative values for concentration.

Is there any way that I can constrain the solutions produced by bvp4c (or bvp5c) to reflect the constraints on the physical entities that they represent? In the present case, C_Ga and C_As represent concentrations so I would like to introduce a non-negativity constraint. I know that the Matlab function "assume" does something like this for symbolic variables, but as far as I can see I cannot make use of this within bvp4c.

Thank you in advance for your help,

Amy

[mosi nezhadhasan]

hi, you can use bvp4c. for example i solve system of couple ode(BVP) :

• y1''-y2=0 with Bc. y1(0)=0 and y1(1)=2
• y2''-1=0 with BC. y2(0)=0 and y2(1)=1

this system of couple BVP has a simple algebraic answer. for this system we can write :

if true
  function SimpleBVP()

     solinit = bvpinit(linspace(0,1,10),[0 0 0 0]);

     sol = bvp4c(@deriv,@bcs,solinit);

     plot(sol.x,sol.y(1,:),'b-x',sol.x,sol.y(2,:),'b-x');

      function dYdx = deriv(x,Y)

       dYdx(1) = Y(3); 
      dYdx(2) = Y(4);
      dYdx(3) = Y(2);
      dYdx(4) = 1;

     % boundary conditions y'(a)=1, y(b)=0.
     function res = bcs(ya,yb)
        res = [ ya(1)
            yb(1)-2
        ya(2)
        yb(2)-1];
  end

[Jasper Wong]

Hi Amy, did you ever resolve this problem? I am facing similar difficulties.

Thanks

[Rodrigo Jorge Leonardi]

Hi Jasper,

¿Did you ever resolve this problem? I am having similar difficulties.

Thanks you.

Answer 1

[John D'Errico]

The simple fix that treats the problem in a way that is arguably correct, is to NOT work directly in terms of concentration! That is, you want to work in a transformed domain. The differential equation solver CANNOT handle constraints on the solution. It does not understand them, and what happens is you get negative numbers. In terms of the mathematics of the differential equation, it makes complete sense. But in terms of what you are modeling, the model fails. So how do we fix this?

Use this transformation into a log domain.

C_Ga = exp(L_Ga)

C_As = exp(L_As)

The unknowns you will solve for are L_Ga and L_As. Just substitute into the problem, for example:

d(C_Ga)/dr = exp(L_GA) * d(L_ga)/dr

The second derivative is as easily generated. Completely replace every term that has C_Ga in it, as well as C_As. Then toww it into a BVP solver.

As you should see, the solution will no longer have any problems with negative concentrations. They no longer exist, since exp(L_Ga) will never be negative.

Once you have the solution in terms of L_Ga and L_As, exponentiate to get the desired concentrations.

[abs mhr]

Dear John What function should be defined to handle the transformation in case of constraints such as -b <C_As< a and -k <C_Ga< l ? Is this solution applicable both for costate and state variables? Thank you

Answer 2

[Jacob Miller]

Hi everyone, I think I may have figured this out. I have a similar issue (also solving a diffusion-reaction boundary value problem, BVP5C was giving me negative concentrations) and came across this post discussing why solvers using a collocation method (like BVP4C/BVP5C) struggle with stiff differential equations, in which a quantity changes rapidly near one boundary, then stays relatively constant throughout the rest of the domain: https://math.stackexchange.com/questions/786148/due-to-numerical-inaccuracy-the-solution-of-a-boundary-value-problems-becomes-n

I changed the mesh I was feeding to BVP5C from one that was linearly spaced to one that was logarithmically spaced. Admittedly, I originally placed most of my points nearest to the outer boundary (where the variables are rapidly changing) in order to give MATLAB more initial guesses at a solution here. This actually didn't work--I instead have most of my points bunched up in the flat area. This got rid of my negative concentration issue. Hope it helps!

[John D'Errico]

See my answer, which explains how to transform the problem so you have no need to play with the meshes, etc.