nonlin

Introduction

NONLIN is a library that provides routines for solving least squares problems, systems of equations, and single equations. The library provides function optimization routines. There are also specific routines and data types optimized for operation on polynomials.

The example below illustrates two techniques of determining the roots of a polynomial using routines and data types provided by this library.

program example
    use iso_fortran_env
    use nonlin_core
    use nonlin_solve
    implicit none
 
    ! Local variables
    type(fcn1var_helper) :: obj
    procedure(fcn1var), pointer :: fcn
    type(newton_1var_solver) :: solver
    real(real64) :: x, f
    type(value_pair) :: limits
    type(iteration_behavior) :: tracking
 
    ! Define the search limits
    limits%x1 = 2.0d0
    limits%x2 = -2.0d0
 
    ! Establish the function
    fcn => fcn1
    call obj%set_fcn(fcn)
 
    ! Allow the solver to print out updates at each iteration
    call solver%set_print_status(.true.)
 
    ! Solve the equation
    call solver%solve(obj, x, limits, f, ib = tracking)
 
    ! Print the output, residual, and details regarding the iteration
    print *, ""
    print '(AF7.5)', "The solution: ", x
    print '(AE10.3)', "The residual: ", f
    print '(AI0)', "Iterations: ", tracking%iter_count
    print '(AI0)', "Function Evaluations: ", tracking%fcn_count
    print '(AI0)', "Derivative Evaluations: ", tracking%jacobian_count
    print '(AL1)', "Converge on Function Value: ", tracking%converge_on_fcn
    print '(AL1)', "Converge on Change in Variable: ", tracking%converge_on_chng
    print '(AL1)', "Converge on Derivative: ", tracking%converge_on_zero_diff
 
contains
    ! The function:
    ! f(x) = x**3 - 2 * x - 1
    function fcn1(x) result(f)
        real(real64), intent(in) :: x
        real(real64) :: f
        f = x**3 - 2.0d0 * x - 1.0d0
    end function
end program

The above program produces the following output.

Iteration: 1
Function Evaluations: 4
Jacobian Evaluations: 2
Change in Variable: 0.500E+00
Residual: -.125E+00
 
Iteration: 2
Function Evaluations: 5
Jacobian Evaluations: 3
Change in Variable: 0.125E+01
Residual: -.208E+01
 
Iteration: 3
Function Evaluations: 6
Jacobian Evaluations: 4
Change in Variable: 0.625E+00
Residual: -.115E+01
 
Iteration: 4
Function Evaluations: 7
Jacobian Evaluations: 5
Change in Variable: -.313E+00
Residual: 0.436E+00
 
Iteration: 5
Function Evaluations: 8
Jacobian Evaluations: 6
Change in Variable: 0.665E-01
Residual: 0.221E-01
 
Iteration: 6
Function Evaluations: 9
Jacobian Evaluations: 7
Change in Variable: 0.375E-02
Residual: 0.685E-04
 
The solution: 1.61803
The residual:  0.665E-09
Iterations: 7
Function Evaluations: 11
Derivative Evaluations: 8
Converge on Function Value: T
Converge on Change in Variable: F
Converge on Derivative: F

An alternative for determining the roots of a polynomial is to use the polynomial type from this library, and determine the roots using the functionallity of that type. The following example illustrates just such an approach.

program example
    use iso_fortran_env
    use nonlin_polynomials
    implicit none
 
    ! Local Variables
    type(polynomial) :: f
    real(real64) :: coeffs(4)
    complex(real64), allocatable :: rts(:)
    integer(int32) :: i
 
    ! Define the polynomial (x**3 - 2 * x - 1)
    coeffs = [-1.0d0, -2.0d0, 0.0d0, 1.0d0]
    f = coeffs
 
    ! Compute the polynomial roots
    rts = f%roots()
 
    ! Display the results
    do i = 1, size(rts)
        print '(AI0AF9.6AF9.6A)', "Root ", i, " = (", real(rts(i), real64), &
            ", ", aimag(rts(i)), ")"
    end do
end program

The above program produces the following output.

Root 1 = (-1.000000,  0.000000)
Root 2 = (-0.618034,  0.000000)
Root 3 = ( 1.618034,  0.000000)