Run "if" logic on a row

I have the following code for the reaction interface between acid and base (A-acid, B-base,S-salt):

    dS=k1*cA(i+1,:).*cB(i+1,:)*dt;
    IDX = (cA(i+1,:)<=NC | cB(i+1,:)<=NC);   % Logical indexing
    %if in a column for this time iteration either cA or cB is less than NC
    %then set the value of dS equal to zero (i.e. no reaction occurred at
    %that location
    dS(IDX)=0;  

where cA and cB are arrays, dS is a vector with the same number of columns as cA and cB and dt is the time increment.

Now I want to analyze each element of the current row (i+1, the time iteration) in cA and cB based on what the corresponding value is in the column (space increment) of the vector dS. If dS equals zero, I want to set cA and cB to one set of values (no reaction has occurred). If dS is greater than zero, then reaction occurred and a different line of code is used to set the value of the concentration in that cell in either cA or cB.

If dS > 0

    cA(i+1,:)=cA(i+1,:)-dS;
    cB(i+1,:)=cB(i+1,:)-dS; 

if dS = 0, seems to be more tricky. In this case the acid and base concentrations are set by the water dissociation constant, Kw. The code needs to determine which array (cA or cB) has the higher value for the concentration. Then the concentration of the one with the lower value is calculated based on Kw and the higher concentration. So if cA = .02 and cB = .0000001 in equivalent cells in cA and cB (equivalent based on row/column), then cA = 0.02 and cB = kW/cA.

Coming from a visual basic (VB) background, I have a hard time understanding how to handle many cells at once as MatLab can do. Clearly, it's more efficient and I don't want to use VB logic to code this issue.

Answer 1

[Joel Lynch]

The trick is to recognize that you can use & / | (and/or) to chain multiple logical operations together. In this case, you can get logical inidices that match both dS=0 and cA>cB.

% indicies of columns where dS==0 AND cA>cB
idx = dS==0 & cA(i+1,:)>cB(i+1,:);
cB(i+1,idx) =  kW/cA(i+1,idx);

% inidices of colums where dS==0 AND cB>cA
idx = dS==0 & cA(i+1,:)<cB(i+1,:);
cA(i+1,idx) =  kW/cB(i+1,idx);

The various operators used in matlab are summarized here, and can be useful when learning the language.

[Robert Demyanovich]

Matlab is throwing up an error at the following line:

cB(i+1,idx) = Kw/cA(i+1,idx);

The error message is:

"Error using /

Matrix dimensions must agree."

cA and cB have the exact same dimensions.

[Joel Lynch]

try element-wise operators "./" instead of "/"

[Robert Demyanovich]

That fixed that error.

But as part of debugging, because the overall results are incorrect, I decided to turn dS into a matrix of the same size as cA and cB so that I could check how the values of dS are changing with the time step.

Now I am getting a different error:

"The logical indices in position 2 contain a true value outside of the array bounds.

Error in TestRJDacidBasePDE (line 213)

cB(i+1,idx) = Kw./cA(i+1,idx);"

This error looks tricky because if I put the cA lines ahead of the cB lines in the code that your provided, I still get the same error at the cB line.

In other words,

idx = dS==0 & cA(i+1,:)<cB(i+1,:);

cA(i+1,idx) = Kw./cB(i+1,idx);

does not generate an error if it precedes:

idx = dS==0 & cA(i+1,:)>cB(i+1,:);

cB(i+1,idx) = Kw./cA(i+1,idx);

Only cB(i+1,idx) = Kw./cA(i+1,idx); generates an error??!

Also, I made the following change in the code I originally posted, changing dS to dS(i+1,:).

    dS(i+1,:)=k1*cA(i+1,:).*cB(i+1,:)*dt;
    IDX = (cA(i+1,:)<=NC | cB(i+1,:)<=NC);   % Logical indexing
    %if in a column for this time iteration either cA or cB is less than NC
    %then set the value of dS equal to zero (i.e. no reaction occurred at
    %that location
    dS(IDX)=0;  

[Joel Lynch]

If dS is now a matrix, you have to include the i+1 row in your logical operations. If you don't, then in the following line

idx = dS==0 & cA(i+1,:)<cB(i+1,:);

should throw an error because dS==0 and cA(i+1,:)<cB(i+1,:) return differently-sized matricies/vectors. The same is true for setting dS(IDX)=0.

So you should have something like:

dS(i+1,:)=k1*cA(i+1,:).*cB(i+1,:)*dt;
IDX = (cA(i+1,:)<=NC | cB(i+1,:)<=NC);   % Logical indexing
%if in a column for this time iteration either cA or cB is less than NC
%then set the value of dS equal to zero (i.e. no reaction occurred at
%that location
dS(i+1,IDX)=0; %<THIS CHANGED
    
% indicies of columns where dS==0 AND cA>cB
idx = dS(i+1,:)==0 & cA(i+1,:)>cB(i+1,:); % < THIS CHANGED
cB(i+1,idx) =  kW./cA(i+1,idx);

% inidices of colums where dS==0 AND cB>cA
idx = dS(i+1,:)==0 & cA(i+1,:)<cB(i+1,:); % < THIS CHANGED
cA(i+1,idx) =  kW./cB(i+1,idx);

Are you still getting the error? If so, can you attach your full script + data?

Out of curisoity, are you time stepping a chemical reaction? If so, are you aware that MATLAB has built-in stiff ODE solvers that are much more stable than euler integration?

[Robert Demyanovich]

I had a feeling that was the problem. I actually did change the first dS line to include i+1, but not the others. So, at least I'm starting to get a feel for how this works.

No more errors, but the results are now way off. I need to figure out why my logic is wrong.

Yes, I am stepping through a chemical reaction that is virtually instantaneous. When I assume that the concentrations at the reaction interface are zero (because of the instantaneous rxn), I get the results that are to be expected. I'm now trying to change that to include Kw, and the results are as I said way off. I'm surprised because I'm trying to set the concentrations at the reaction plane to 1E-07 (pH =7) for both acid and base instead of concentrations of zero. Since the starting concentrations are 5 orders of magnitude higher, I would have thought this is pretty darn close to zero so the results should be almost identical. Must be a logic error somewhere.

I was not aware of the built-in stiff ODE solver. It's true that if I am not careful about the number of steps chosen for both time and space, the algorithm blows up. I guess the ODE solver would fix this. I only purchased the bare bones MatLab - so I'm not sure if the built in solver is available to me.

[Joel Lynch]

ODE solvers should be included in the base software. MATLAB has some good documentation on how to choose and implement an ode scheme. I would be vary careful using euler integration, unless you know the system is non-stiff.

Questions on chemical kinetics are probably outside the scope of this forum, but a standard way to test stiff ODE's is to continually reduce your time-step size. If your value converges to a constant result, stiffness probably isn't the most immediate problem.

Good luck!