# Newton-Raphson routine

iherrera

Mon, 10/03/2011 - 02:09 pm

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

// 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 />

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

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

Wave w, xW, yW

yW[0] = w[0][0]*xW[0]+w[0][1]*xW[0]*xW[1]+2*w[0][2]*xW[0]^2*xW[1] -w[0][4]

yW[1] = w[1][0]*xW[1]+w[1][1]*xW[0]*xW[1]+ w[1][2]*xW[0]^2*xW[1]+ w[1][3]*xW[1]*xW[2]-w[1][4]

yW[2] = w[2][0]*xW[2]+ w[2][3]*xW[1]*xW[2]-w[2][4]

End

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