Solver, Minimization, 2 Variables

[ikkebins]

Hey there,

following problem:

A set of non linear equations, that can be solved by determining two variables.

So far i have some experience using the NEQNF Solver, that minimizes an objective function (or several functions) by varying one variable.

Now i have 7 objective functions (non linear) and i need to vary and determine two variables to solve minimize the functions. Is there an solver that you can suggest me to use for this problem?


I`m using fortrann 77.

Thanks for your help,
ikkebins

 

[salgerman]

I don't know; but, have you looked for Optimization in netlib?

 

[ikkebins]

Yes, but i couldnt find anything

 

[mikrom]

If you cannot find nothing appropriate, you can try to code it yourself.
For solving the nonlinear equation system of this form

f(x) = 0

we can use Newton method:

x[sub]1[/sub] = x[sub]0[/sub] - J[sup]-1[/sup](f, x[sub]0[/sub]) * f(x[sub]0[/sub])

That means:

1) For the given starting point x[sub]0[/sub] we need to do:

1.a) Compute the Jacobian matrix J(f, x[sub]0[/sub]) of the vector function f(x) in the point x[sub]0[/sub]

1.b) Solve the linear equation system
J(f, x[sub]0[/sub]) * d = - f(x[sub]0[/sub])
For this step we can use a linear equation solver from Lapack.

2) finaly compute the next approximation point
x[sub]1[/sub] = x[sub]0[/sub] + d

All steps should be repeated until the desired accuracy is reached, it means for given eps -> 0 or delta -> 0 until
||x[sub]n[/sub] - x[sub]n-1[/sub]|| < delta or ||f(x[sub]n[/sub])|| < eps

IMO it doesn't seem very complicated...

 

[mikrom]

I looked at the description of NEQNF method here

http://faculty.ksu.edu.sa/Almutaz/Documents/Summer-2010/ChE-401/NEQNF.pdf

I took as example the 3x3 system of nonlinear equations from the document above and tried to solve it with simple Newton method.
For solving of linear equation systems I first used self written GEM method, but then - to make the code shorter for this posting - I took a method from LAPACk:

nlsolve.f95

[COLOR=#a020f0]module[/color] linear_algebra
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
  [COLOR=#0000ff]! This module contains the tools[/color]
  [COLOR=#0000ff]! for Linear Algebra[/color]
[COLOR=#a020f0]contains[/color]
[COLOR=#2e8b57][b]  real[/b][/color] [COLOR=#a020f0]function[/color] norm(X)
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](:) :: X
    [COLOR=#0000ff]! short notation of Fortran 90[/color]
    norm [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]sqrt[/color]([COLOR=#008080]sum[/color](X[COLOR=#804040][b]*[/b][/color]X))
  [COLOR=#a020f0]end function[/color] norm

  [COLOR=#a020f0]subroutine[/color] Lsolve(M, V, X, n)
    [COLOR=#0000ff]! LAPACK routine for solving linear equation system[/color]
    [COLOR=#0000ff]! M * X = V[/color]
    [COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: n
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n,n), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]inout[/b][/color]) :: M
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]inout[/b][/color]) :: V
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]out[/b][/color]) :: X 
    [COLOR=#0000ff]![/color]
[COLOR=#2e8b57][b]    real[/b][/color] :: A(n,n), b(n)
    [COLOR=#2e8b57][b]integer[/b][/color] :: i, pivot(n), ok
    A [COLOR=#804040][b]=[/b][/color] M
    b [COLOR=#804040][b]=[/b][/color] V
    [COLOR=#008080]call[/color] SGESV(n, [COLOR=#ff00ff]1[/color], A, n, pivot, b, [COLOR=#ff00ff]3[/color], ok)    
    X [COLOR=#804040][b]=[/b][/color] b
  [COLOR=#a020f0]end subroutine[/color] Lsolve
[COLOR=#a020f0]end module[/color] linear_algebra

[COLOR=#a020f0]module[/color] solving_methods
  [COLOR=#a020f0]use[/color] linear_algebra
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
  [COLOR=#0000ff]! This module contains the solving method[/color]
[COLOR=#a020f0]contains[/color]
  [COLOR=#a020f0]subroutine[/color] jacobian(F, P, Jcb, n)
    [COLOR=#2e8b57][b]integer[/b][/color] :: n
    [COLOR=#0000ff]! interface for functional argument[/color]
    [COLOR=#a020f0]interface[/color]
      [COLOR=#a020f0]function[/color] F(X, n)
        [COLOR=#2e8b57][b]integer[/b][/color] :: n
[COLOR=#2e8b57][b]        real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: x
[COLOR=#2e8b57][b]        real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n) :: f
      [COLOR=#a020f0]end function[/color] f
    [COLOR=#a020f0]end interface[/color]    
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: P
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n,n), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]out[/b][/color]) :: Jcb
    [COLOR=#0000ff]![/color]
    [COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]parameter[/b][/color] :: info [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color] [COLOR=#0000ff]! view informations[/color]
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]parameter[/b][/color] :: h[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]0.1[/color]   [COLOR=#0000ff]! step    [/color]
    [COLOR=#2e8b57][b]integer[/b][/color] :: i, j
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n) :: diffj, T1, T2
    [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], n   [COLOR=#0000ff]! nr of functions[/color]
      [COLOR=#804040][b]do[/b][/color] j[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], n [COLOR=#0000ff]! nr of variables[/color]
        T1 [COLOR=#804040][b]=[/b][/color] P
        T2 [COLOR=#804040][b]=[/b][/color] P
        T1(j) [COLOR=#804040][b]=[/b][/color] P(j) [COLOR=#804040][b]+[/b][/color] h
        T2(j) [COLOR=#804040][b]=[/b][/color] P(j) [COLOR=#804040][b]-[/b][/color] h

        diffj [COLOR=#804040][b]=[/b][/color] (F(T1, n) [COLOR=#804040][b]-[/b][/color] F(T2, n)) [COLOR=#804040][b]/[/b][/color] ([COLOR=#ff00ff]2[/color][COLOR=#804040][b]*[/b][/color]h)
        [COLOR=#804040][b]if[/b][/color] (info [COLOR=#804040][b].ne.[/b][/color] [COLOR=#ff00ff]0[/color]) [COLOR=#804040][b]then[/b][/color]
           [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'T1='[/color], T1
           [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'T2='[/color], T2
           [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'diffj = '[/color], diffj
        [COLOR=#804040][b]end if[/b][/color]
        [COLOR=#0000ff]! difference of i-th function on j-th variable[/color]
        Jcb(i,j) [COLOR=#804040][b]=[/b][/color] diffj(i)
      [COLOR=#804040][b]end do[/b][/color]
    [COLOR=#804040][b]end do[/b][/color]
  [COLOR=#a020f0]end subroutine[/color] jacobian

  [COLOR=#a020f0]subroutine[/color] newton(F, X, n)
    [COLOR=#0000ff]! Newton-Raphson Method for solving nonlinear systems[/color]
    [COLOR=#2e8b57][b]integer[/b][/color] :: n
    [COLOR=#0000ff]! interface for functional argument[/color]
    [COLOR=#a020f0]interface[/color]
      [COLOR=#a020f0]function[/color] F(X, n)
        [COLOR=#2e8b57][b]integer[/b][/color] :: n
[COLOR=#2e8b57][b]        real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: x
[COLOR=#2e8b57][b]        real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n) :: f
      [COLOR=#a020f0]end function[/color] f
    [COLOR=#a020f0]end interface[/color]    
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]inout[/b][/color]) :: X
    [COLOR=#0000ff]![/color]
    [COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]parameter[/b][/color] :: info [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1[/color] [COLOR=#0000ff]! view informations[/color]
    [COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]parameter[/b][/color] :: max_steps [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]100[/color]
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]parameter[/b][/color] :: delta[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1e-5[/color], eps[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1e-8[/color]
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n,n) :: JF
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n) :: T, D, X_old
[COLOR=#2e8b57][b]    real[/b][/color] :: dnorm, fnorm
    [COLOR=#2e8b57][b]integer[/b][/color] :: k
    
    [COLOR=#804040][b]do[/b][/color] k[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], max_steps
      X_old [COLOR=#804040][b]=[/b][/color] X
      T [COLOR=#804040][b]=[/b][/color] F(X, n)
      [COLOR=#008080]call[/color] Jacobian(F,X,JF,n)
      [COLOR=#008080]call[/color] Lsolve(JF,T,D,n)
      X [COLOR=#804040][b]=[/b][/color] X [COLOR=#804040][b]-[/b][/color] D
      dnorm [COLOR=#804040][b]=[/b][/color] norm(D)
      fnorm [COLOR=#804040][b]=[/b][/color] norm(T)
      [COLOR=#804040][b]if[/b][/color] ((fnorm [COLOR=#804040][b]<=[/b][/color] eps) [COLOR=#804040][b].or.[/b][/color] (dnorm [COLOR=#804040][b]<=[/b][/color] delta)) [COLOR=#804040][b]then[/b][/color]
        [COLOR=#0000ff]! exit loop [/color]
        [COLOR=#804040][b]exit[/b][/color]
      [COLOR=#804040][b]end if[/b][/color]      
    [COLOR=#804040][b]end do[/b][/color]
    [COLOR=#804040][b]if[/b][/color] (info [COLOR=#804040][b].ne.[/b][/color] [COLOR=#ff00ff]0[/color]) [COLOR=#804040][b]then[/b][/color]
      [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'fnorm ='[/color], fnorm
      [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'dnorm ='[/color], dnorm
      [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'nr steps ='[/color], k
    [COLOR=#804040][b]end if[/b][/color]    
  [COLOR=#a020f0]end subroutine[/color] newton
[COLOR=#a020f0]end module[/color] solving_methods

[COLOR=#a020f0]module[/color] user_functions
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
  [COLOR=#0000ff]! Define here your own right-hand function[/color]
  [COLOR=#0000ff]! for the non linear system[/color]
  [COLOR=#0000ff]!    F(X) = 0[/color]
[COLOR=#a020f0]contains[/color]
  [COLOR=#a020f0]function[/color] nls1(X, n)
    [COLOR=#0000ff]! starting point: 4.0, 4.0, 4.0[/color]
    [COLOR=#0000ff]! solution      : 1.0 2.0 3.0[/color]
    [COLOR=#2e8b57][b]integer[/b][/color] :: n
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: X
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n) :: nls1
    [COLOR=#0000ff]![/color]
    NLS1([COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]=[/b][/color] X([COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]+[/b][/color] [COLOR=#008080]EXP[/color](X([COLOR=#ff00ff]1[/color])[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1.0[/color]) [COLOR=#804040][b]+[/b][/color] (X([COLOR=#ff00ff]2[/color])[COLOR=#804040][b]+[/b][/color]X([COLOR=#ff00ff]3[/color]))[COLOR=#804040][b]*[/b][/color](X([COLOR=#ff00ff]2[/color])[COLOR=#804040][b]+[/b][/color]X([COLOR=#ff00ff]3[/color])) [COLOR=#804040][b]-[/b][/color] [COLOR=#ff00ff]27.0[/color]
    NLS1([COLOR=#ff00ff]2[/color]) [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]EXP[/color](X([COLOR=#ff00ff]2[/color])[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]2.0[/color])[COLOR=#804040][b]/[/b][/color]X([COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]+[/b][/color] X([COLOR=#ff00ff]3[/color])[COLOR=#804040][b]*[/b][/color]X([COLOR=#ff00ff]3[/color]) [COLOR=#804040][b]-[/b][/color] [COLOR=#ff00ff]10.0[/color]
    NLS1([COLOR=#ff00ff]3[/color]) [COLOR=#804040][b]=[/b][/color] X([COLOR=#ff00ff]3[/color]) [COLOR=#804040][b]+[/b][/color] [COLOR=#008080]SIN[/color](X([COLOR=#ff00ff]2[/color])[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]2.0[/color]) [COLOR=#804040][b]+[/b][/color] X([COLOR=#ff00ff]2[/color])[COLOR=#804040][b]*[/b][/color]X([COLOR=#ff00ff]2[/color]) [COLOR=#804040][b]-[/b][/color] [COLOR=#ff00ff]7.0[/color]
  [COLOR=#a020f0]end function[/color] nls1
[COLOR=#a020f0]end module[/color] user_functions

[COLOR=#a020f0]program[/color] nlsys
  [COLOR=#a020f0]use[/color] solving_methods
  [COLOR=#a020f0]use[/color] user_functions
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
  
[COLOR=#2e8b57][b]  real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color]([COLOR=#ff00ff]3[/color]) :: X

  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Solving Nonlinear System:'[/color]
  x [COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]4[/color],[COLOR=#ff00ff]4[/color],[COLOR=#ff00ff]4[/color][COLOR=#804040][b]/[/b][/color])
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Starting Point ='[/color], X
  [COLOR=#008080]call[/color] newton(nls1, X, [COLOR=#ff00ff]3[/color])
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Solution ='[/color], X
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])
[COLOR=#a020f0]end program[/color] nlsys

Here are the results:

$ gfortran -o nlsolve nlsolve.f95 -L. liblapack.dll

$ nlsolve
 Solving Nonlinear System:
 Starting Point =   4.00000000       4.00000000       4.00000000
 fnorm =   1.33280039E-07
 dnorm =   5.93307874E-08
 nr steps =           7
 Solution =   1.00000000       2.00000000       3.00000000