Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
fsolve for competing equilibria

Subject: fsolve for competing equilibria

From: Jesper

Date: 13 Mar, 2012 16:46:13

Message: 1 of 2

I'm trying to solve an equilibrium problem for the following competing chemical reaction:
AP + P -><- APa@P (K1 = APa@P/(AP*P) = 1e5)
AP + P -><- APp@P (K2 = APp@P/(AP*P) = 1e4)

That is; AP can bind in two ways to P and does so with two different association constants (K1 and K2). I will therefore have 4 species in solution: unbound AP, unbound P, APa@P and APp@P. Knowing the starting concentrations AP0 and P0, fsolve should easily be able to solve all the concentrations by solving the following four non-linear equations:

P0 - P - APa@P - APp@P = 0
AP0 - APa@P - APp@P = 0
APa@P/(P*AP) - K1 = 0
APp@P/(P*AP) - K2 = 0

My function file becomes:

function [Concentrations] = svdConcfunc(D)

global AP0 P0 K1 K2

   Concentrations = [P0 - D(2) - D(3) - D(4);
                    AP0 - D(1) - D(3) - D(4);
                    D(3)./(D(1).*D(2)) - K1;
                    D(4)./(D(1).*D(2)) - K2];
 return

My script, calling the above function looks like this:

clear all
close all
format long

global conc AP0 P0 K1 K2

options = optimset ('Tolfun', 1e-30, 'TolX', 1e-20, 'MaxFunEvals', 1e20, 'MaxIter', 1e20);

conc=load('AP_titration_conc_vector.txt')./(1e6);

P0 = 2e-6;
K1 = 1e5;
K2 = 1e4;

for i = 1:length(conc);
AP0 = conc(i);
%format = [AP P APa@P APp@P ];
Guessconc = [AP0 P0 AP0 AP0 ];
CalcConc(i,:) = fsolve(@svdConcfunc, Guessconc, options);
test_should_be(i) = P0 + AP0;
end

test = sum(CalcConc, 2)
test_should_be = test_should_be'

the loaded "conc" becomes:
conc=

  1.0e-004 *

                   0
   0.000800000000000
   0.001600000000000
   0.006380000000000
   0.014300000000000
   0.030000000000000
   0.053000000000000
   0.082800000000000
   0.126000000000000
   0.180000000000000
   0.244000000000000
   0.315000000000000
   0.396000000000000
   0.493000000000000

The output I get is the following:

test =

  1.0e-004 *

   0.020000000000000
   0.020656586717874
   0.021314855866978
   0.025281807267153
   0.031972557259257
   0.045604844303509
   0.066255477286183
   0.093836312858125
   0.134837016836453
   0.187048378261830
   0.249671527576348
   0.319655268142224
   0.399862134556760
   0.496206019001146


test_should_be =

  1.0e-004 *

   0.020000000000000
   0.020800000000000
   0.021600000000000
   0.026380000000000
   0.034300000000000
   0.050000000000000
   0.073000000000000
   0.102800000000000
   0.146000000000000
   0.200000000000000
   0.264000000000000
   0.335000000000000
   0.416000000000000
   0.513000000000000

So the sum of all calculated concentrations in solution is not the same (there is a difference between "test_should_be" and "test") as the sum of what I put in, which is physically impossible and I cannot for the life of me understand why it won't work. I have been able to solve far more complicated equilibrium problems (more equations, species and binding constants) before so this should be a walk in the park.

I have tried to play around with starting guesses and tolerance values but I can't seem to get a better solution.

Any help on this would be much appreciated!

/Jesper

Subject: fsolve for competing equilibria

From: Roger Stafford

Date: 13 Mar, 2012 21:06:19

Message: 2 of 2

"Jesper" wrote in message <jjntkl$8p1$1@newscl01ah.mathworks.com>...
> .......
> test_should_be(i) = P0 + AP0;
> .......
> test = sum(CalcConc, 2)
> ........
> So the sum of all calculated concentrations in solution is not the same (there is a difference between "test_should_be" and "test") as the sum of what I put in, which is physically impossible and I cannot for the life of me understand why it won't work. ....
- - - - - - - - -
  Though I am not a chemist there seems to be something rather questionable with your test procedure. For every one mole of AP and one mole of P which is converted to APa@P or APp@P you have only one mole of result, but you appear to be requiring in "test=sum(CalcConc,2)" that the total number of moles should remain constant. In other words you may be obtaining perfectly valid solutions but your test looks invalid to me.

  A better test would be to check whether your original four equations are satisfied or not for each of the AP0 concentrations.

  There is a way you can avoid using 'fsolve' and I don't understand why you are not doing things this way. Your four equations can be reduced to the solution of simple quadratic equations in D1 and D2 (AP and P) in which one of the two roots gives positive values for the concentrations. D3 and D4 (APa@P and APp@P) are then readily found from D1 and D2. Do as follows:

 T1 = (-1+sqrt(1+2*(K1+K2)*(P0+AP0)+((K1+K2)*(P0-AP0)).^2))/2/(K1+K2);
 T2 = (P0-AP0)/2;
 D1 = T1-T2;
 D2 = T1+T2;
 D3 = K1*D1.*D2;
 D4 = K2*D1.*D2;

(These should be valid as vectorized solutions using a vector for AP0.)

  I am assuming I have the correct equations to be solved:

 D1 + D3 + D4 = AP0
 D2 + D3 + D4 = P0
 D3/(D1*D2) = K1
 D4/(D1*D2) = K2

(You first wrote: "AP0-APa@P-APp@P=0" for that first equation but later added the obviously needed D1 (AP) term.)

Roger Stafford

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us