Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

Register Log in
  • Congratulations SkipVought on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

Solver, Minimization, 2 Variables

Status
Not open for further replies.
ikkebins

ikkebins

Technical User
May 28, 2014
2
0
0
DE
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
 
Sort by date Sort by votes
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...
 
Upvote 0 Downvote
I looked at the description of NEQNF method here
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
Code: 
[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:
Code: 
$ 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
 
Upvote 0 Downvote
Status
Not open for further replies.

Similar threads

ksqr
  • Locked
  • Question
Need a 2D vector graphics package
Replies
5
Views
150
mikrom
mikrom
WAX13
  • Locked
  • Question
file extension for free form FORTRAN 2018
Replies
4
Views
87
mikrom
mikrom
JManuelR
  • Locked
  • Question
Sharing variables as an argument with Modules 1
Replies
7
Views
58
JManuelR
JManuelR
COQ
  • Locked
  • Question
Compiling MVS 2010 Fortran project with recourse to .mod and .o module files
Replies
1
Views
59
xwb
xwb
mikkelstuff
  • Locked
  • Question
gfortran: loading hex into DATA statement CHARACTER variables
Replies
8
Views
53
xwb
xwb

Part and Inventory Search

Sponsor

Back
Top