Newton-Raphson routine
iherrera
I'm trying to implement a Newton-Raphson routine to solve chemical equilibrium . I found some code for this routine in Matlab and I've translated it into Igor, but at the moment my routine fails to minimize. I've tested the code and everything works well except for the calculation of the shift vector. The shift vector is computed with the lines
MatrixOp /O mDelta = Inv(mJcb)*Diagonal(Cf) and delta_C = dif_C[p]*mDelta[p][p]*FxC . I'm not sure if I've translated this code properly. I would appreciated if anyone with experience in Matlab and Igor could help me sort this out. For reference, I've posted my code and the code written for Matlab.
I would also appreciated if anyone has suggestions on how to solve this same problem with the FindRoots function in Igor. The mass balance equations that I'm trying to solve have the general form for components [A],[B], and [H] as
At= [A] + k110[A][B] + 2*k210[A]^2[B]
Bt= [B] + k110[A][B] + k210[A]^2[B]+k011[B][H]
Ht= [H] + k011[B][H]
or in matrix form m=
1 0 0 1 2 0
0 1 0 1 1 1
0 0 1 0 0 1
k = {1,1,1,k110,k120,k011}
Ct= {At,Bt,Ht}
fxC={1,1,0} //Constant pH
Here's the code I wrote in Igor
Function NewtonRaphson(m, k , Ct , Ci, fxC, It) // m = model in matrix form // k = binding constants for species (i.e. K=[HX]/([H][X]) , by definitition K=1 for the free components, (i.e. K= [H]/[H]=1) // Ct = total concentration for Components // Ci = Initial guess values for Components // FxC =wave used to fix concentration of components in Co (i.e. fx={1,1,1,0,1}, fourth component's concentration is fixed) wave m, k , Ct , Ci, fxC variable It Variable n_comp = numpnts (Ct) //number of componets Variable n_spec = numpnts (k) //numer of species Variable ii = 0 Variable jj,kk=0 //Check for zero values in Co and replace them with 1e-15 Ct=Ct[p]==0?1e-15:Ct[p] Duplicate /Free k C_spec, C_Jcb Duplicate /Free Ct calc_Ct, dif_C, sq_dif, delta_C,absD_C Duplicate /Free Ci Cf Duplicate /Free m wM //Dummy matrix of similar dimentions to model used to perform matrix calculations Make /Free /N=(n_comp,n_comp), mJcb, mDelta //square Jacobian and delta matrix do ii+=1 //Concentration for species wM = Cf[p] wM = wM^m C_spec = Prod(wM,1)*K //Calculated total concentration for components wM = C_spec[q] wM = m*wM MatrixOp /O calc_Ct = sumRows(wM) //Calculate difference between total concentration and calculated concentration for components dif_C = Ct-calc_Ct absD_C = abs(dif_C) if (all(absD_C, "<", 1e-15)) print "small dif" return Cf else //Calculate Jacobian from model matrix and the concentration of species for(jj=0;jj<n_comp;jj+=1) for (kk=jj;kk<n_comp;kk+=1) C_Jcb=m[jj][p]*m[kk][p]*C_spec mJcb[jj][kk]=sum(C_Jcb) mJcb[kk][jj]=mJcb[jj][kk] endfor endfor endif //Calculate shift in concentrations //at the moment the next two lines are not calculating the shift wave properly MatrixOp /O mDelta = Inv(mJcb)*Diagonal(Cf) delta_C = dif_C[p]*mDelta[p][p]*FxC Cf += delta_C if (any(Cf,"<=",0)) do delta_C *= 0.5 Cf = Cf-delta_C absD_C = abs(delta_C) if (all(absD_C, "<", 1e-15)) break endif while (any(Cf,"<=",0)) endif while (ii < It) return Cf End
And here's the code I found for Matlab:
<br /> function c_spec=NewtonRaphson(Model, beta, c_tot, c,i)<br /> <br /> ncomp=length(c_tot); % number of components<br /> nspec=length(beta); % number of species<br /> c_tot(c_tot==0)=1e-15; % numerical difficulties if c_tot=0<br /> <br /> it=0;<br /> while it<=99<br /> it=it+1;<br /> c_spec =beta.*prod(repmat(c',1,nspec).^Model,1); %species conc<br /> c_tot_calc=sum(Model.*repmat(c_spec,ncomp,1),2)'; %comp ctot calc<br /> d =c_tot-c_tot_calc; % diff actual and calc total conc<br /> <br /> if all(abs(d) <1e-15) % return if all diff small<br /> return<br /> end<br /> for j=1:ncomp % Jacobian (J*)<br /> for k=j:ncomp<br /> J_s(j,k)=sum(Model(j,:).*Model(k,:).*c_spec);<br /> J_s(k,j)=J_s(j,k); % J_s is symmetric<br /> end<br /> end<br /> delta_c=(d/J_s)*diag(c); % equation (2.43)<br /> c=c+delta_c;<br /> <br /> while any(c <= 0) % take shift back if conc neg.<br /> delta_c=0.5*delta_c;<br /> c=c-delta_c;<br /> if all(abs(delta_c)<1e-15)<br /> break<br /> end<br /> end<br /> end<br />
I also wrote a couple of subroutines for which I couldn't find an equivalent in Igor. They have limited functionality at the moment... I only wrote code for the cases that I needed.
Function all (w, cond , n) wave w string cond variable n variable np = numpnts (w) variable ii = 0 variable boolean=0 strswitch(cond) // string switch case "<": do if (w[ii] < n) boolean = 1 else boolean = 0 endif ii +=1 while (ii < np && boolean==1) return boolean break case ">": do if (w[ii] > n) boolean = 1 else boolean = 0 endif ii +=1 while (ii < np && boolean==1) return boolean break default: return boolean endswitch return boolean End Function any ( w,cond,n) wave w string cond variable n variable np = numpnts (w) variable ii = 0 variable boolean = 0 strswitch(cond) // string switch case "<=": do if (w[ii] <= n) boolean=1 endif ii +=1 while (ii < np && boolean==0) return boolean break case ">=": do if (w[ii] >= n) boolean = 1 endif ii +=1 while (ii < np && boolean==0) return boolean break default: return boolean endswitch return boolean End Function Prod (w,dim) wave w variable dim variable n = DimSize (w,dim) variable ii, jj Make/O/N=(n)/Free w1=1 if (dim==0) variable col= DimSize (w,1) for(ii=0; ii<n;ii+=1) for(jj=0; jj<col;jj+=1) w1[ii] *= w[ii][jj] endfor endfor elseif (dim==1) variable row= DimSize (w,0) for(ii=0; ii<n;ii+=1) for(jj=0; jj<row;jj+=1) w1[ii] *= w[jj][ii] endfor endfor endif return w1 End
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
October 5, 2011 at 10:18 am - Permalink
and trasfer a matrix as the wave parameter of this form
w^t=
1, k110, k210, 0, At
1, k110, k210, k011, Bt
1, 0, 0, k011, Ht
I tried using FindRoots on the previous function but the values I obtained for the roots have no physical meaning.
IH
October 6, 2011 at 07:55 am - Permalink
Disclaimer: I haven't looked at your code. However, when I look at the equations that you mention in your original post:
I'm assuming that these are rate equations, and that "At" really means "the derivative of the concentration of species A with respect to time".
If that is the case then your equations seem to describe a system in which mass is being created out of nowhere. Typically each growth term, e.g. k110[A][B], should show up with a positive sign in one equation and a negative sign in another one. However, that is not the case in your equations.
I may be misunderstanding your intentions. But it looks odd to me.
October 6, 2011 at 08:59 am - Permalink
Hm... Are you finding roots of a system of equations, or trying to optimize a function? Or maybe you talked about optimizing, because you were trying to roll your own root finder?
In either case, actually, the pWave parameter is a wave that Igor passes along to your function intact without examination. It can be anything you want it to be as long as you interpret it correctly in your function.
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
October 6, 2011 at 10:40 am - Permalink
These equations are mass balance equations for a titration experiment. This means that A(total) represents the sum of the concentrations for the species that contain component A in their formula. For example, the mass balance equations for a system where component A forms complexes with B of type [AB] and [A2B] are defined as:
At= [A]+[AB]+2*[A2B]
Bt= [B]+[AB]+[A2B]
Then, I can express the concentration of [AB] and [A2B] in terms of the equilibrium constants (K) and the concentration of the free components [A] and [B] as:
[AB]= K110[A][B]
[A2B]=K210[A]^2[B]
which leads to the mass balance equations that I wrote before. In that particular system B also coordinates with H, so I added an extra term in Bt and also an additional mass balance equation for Ht.
October 6, 2011 at 10:57 am - Permalink
I am looking for the roots to this system of equations
[A] + k110[A][B] + 2*k210[A]^2[B]- At=0
[B] + k110[A][B] + k210[A]^2[B]+k011[B][H]-Bt=0
[H] + k011[B][H] - Ht =0
which would be the same as writing the polynomials
x + k110xy+2*k210x^2y- At =0
y + k110xy + k210x^2y+k011yz-Bt = 0
z + k011yz-Ht=0
where x=[A], y=[B], and z=[H]. I tried to create my own root finder function because it is very practical to write the chemical equilibrium in matrix form together with a vector that contains the binding constants {1,1,1,k110,k210,k011} and a second vector that contains the total concentration of the components {At,Bt,Ht} (independent variables of my experiments).
At the moment, the roots that I've obtained with FindRoots have no physical meaning because [A] is much larger than At.
October 6, 2011 at 03:10 pm - Permalink
That sounds like a mathematical problem- your system is poorly defined because one component is small.
FindRoots is based on code written by smarter people that me- it has lots of attention paid to things like trying to do the arithmetic in ways that minimize floating-point truncation. Perhaps your system is nearly singular- that would mean you simply can't solve it as it stands, no matter whose code you use. It might be that you need to use one system at the start when the products are at very low concentration, and another at the end where the products are significant.
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
October 7, 2011 at 09:17 am - Permalink
I believe there is nothing wrong with the solver.
There is one missing term in Bt (total concentration of B) that is [BH]
Bt= [B]+[BH]+[AB]+[A2B]
In addition, equation 4 doesn't makes much sense ; [A2B]=K210[A]^2[B]
because it is equivalent to three body collision, but lets assume it holds for the moment.
It should be [A2B]=K210[A][AB], a two body reaction between AB and A
Prof. V. Calvo-Perez
Universidad de Chile
February 19, 2012 at 03:49 pm - Permalink