Newton-Raphson routine

Share icon

Hello everyone,

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
thanks John, I'll give it a try. Would it be OK to write the multivariate function like this

Function ChemEq(w, xW, yW)
  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

iherrera wrote:

I tried using FindRoots on the previous function but the values I obtained for the roots have no physical meaning.


Disclaimer: I haven't looked at your code. However, when I look at the equations that you mention in your original post:
iherrera wrote:

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]

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

iherrera wrote:
and trasfer a matrix as the wave parameter of this form

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

741 wrote:

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".


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

johnweeks wrote:

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?


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

iherrera wrote:
At the moment, the roots that I've obtained with FindRoots have no physical meaning because [A] is much larger than At.

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

iherrera wrote:
741 wrote:

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".


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.



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