- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hello,
Currently, I am using the dtrnlsp_solve (trust region )to sovle a systems of 4 nonlinear equations. I am following the example provided in the https://software.intel.com/en-us/node/471544 Howerver , the powell method failured to get solution.
When I resorted to the fsolve function in Matlab, which uses the trust region method by default as well (pls refer to the link https://www.mathworks.com/help/optim/ug/fsolve.html, I can sucessful get the root of the equations.
Can you pls offer me some suggestions on how to solve it using Dtrnlsp_solver? Thank you in advance。
the matlab code is shown below:
function F = EQs(x) k=190000/3/(1-2*0.3); g=190000/(1+0.3)/2; pe=-167.6485 ; qe=500.0385 ; p=pe+k*x(1); q=qe-3*g*x(2); peeq=x(3); sm=500+1000*(peeq); freal=x(4); ds=sinh(-1.5*p/sm); dc=cosh(-1.5*p/sm); a=0.04*exp(-0.5*((peeq-0.3)/0.1)^2)/((2*3.14159)^0.5)/0.1; F(1) = x(1)* 2*q/sm/sm+x(2)*(-3*1.5*freal*ds/sm); F(2) = (q/sm)^2+2*freal*1.5*dc-(1+2.25*freal^2); F(3)=x(3)-(-p*x(1)+q*x(2))/(1-freal)/sm F(4)=x(4)-(1-freal)*x(1)-a*peeq >>options = optimoptions('fsolve','Display','iter','PlotFcn',@optimplotfirstorderopt); >>x0=0; >>[x,fval] = fsolve(@EQs,x0,options)
the fortran code is shown below :
!=============================================================================== ! Copyright 2004-2016 Intel Corporation All Rights Reserved. ! ! The source code, information and material ("Material") contained herein is ! owned by Intel Corporation or its suppliers or licensors, and title to such ! Material remains with Intel Corporation or its suppliers or licensors. The ! Material contains proprietary information of Intel or its suppliers and ! licensors. The Material is protected by worldwide copyright laws and treaty ! provisions. No part of the Material may be used, copied, reproduced, ! modified, published, uploaded, posted, transmitted, distributed or disclosed ! in any way without Intel's prior express written permission. No license under ! any patent, copyright or other intellectual property rights in the Material ! is granted to or conferred upon you, either expressly, by implication, ! inducement, estoppel or otherwise. Any license under such intellectual ! property rights must be express and approved by Intel in writing. ! ! Unless otherwise agreed by Intel in writing, you may not remove or alter this ! notice or any other notice embedded in Materials by Intel or Intel's ! suppliers or licensors in any way. !=============================================================================== * Content : TR Solver Fortran-77 example * ******************************************************************************** !** NONLINEAR LEAST SQUARE PROBLEM WITHOUT BOUNDARY CONSTRAINTS PROGRAM EXAMPLE_DTRNLSP_POWELL IMPLICIT NONE !** HEADER-FILE WITH DEFINITIONS (CONSTANTS, EXTERNALS) INCLUDE "mkl_rci.fi" !** USER'S OBJECTIVE FUNCTION EXTERNAL EXTENDED_POWELL !** N - NUMBER OF FUNCTION VARIABLES INTEGER N PARAMETER (N = 4) !** M - DIMENSION OF FUNCTION VALUE INTEGER M PARAMETER (M = 4) DOUBLE PRECISION pe,qe,k,g common pe,qe,k,g !** SOLUTION VECTOR. CONTAINS VALUES X FOR F(X) DOUBLE PRECISION X (N) !** PRECISIONS FOR STOP-CRITERIA (SEE MANUAL FOR MORE DETAILS) DOUBLE PRECISION EPS (6) !** JACOBI CALCULATION PRECISION DOUBLE PRECISION JAC_EPS !** REVERSE COMMUNICATION INTERFACE PARAMETER INTEGER RCI_REQUEST !** FUNCTION (F(X)) VALUE VECTOR DOUBLE PRECISION FVEC (M) !** JACOBI MATRIX DOUBLE PRECISION FJAC (M, N) !** NUMBER OF ITERATIONS INTEGER ITER !** NUMBER OF STOP-CRITERION INTEGER ST_CR !** CONTROLS OF RCI CYCLE INTEGER SUCCESSFUL !** MAXIMUM NUMBER OF ITERATIONS INTEGER ITER1 !** MAXIMUM NUMBER OF ITERATIONS OF CALCULATION OF TRIAL-STEP INTEGER ITER2 !** INITIAL STEP BOUND DOUBLE PRECISION RS !** INITIAL AND FINAL RESIDUALS DOUBLE PRECISION R1, R2 !** TR SOLVER HANDLE INTEGER*8 HANDLE !** CYCLE'S COUNTERS INTEGER I, J !** RESULTS OF INPUT PARAMETER CHECKING INTEGER INFO(6) !** SET PRECISIONS FOR STOP-CRITERIA DO I = 1, 6 EPS (I) = 1.D-5 END DO !** SET MAXIMUM NUMBER OF ITERATIONS ITER1 = 1000 !** SET MAXIMUM NUMBER OF ITERATIONS OF CALCULATION OF TRIAL-STEP ITER2 = 100 !** SET INITIAL STEP BOUND RS = 100.D0 !** PRECISIONS FOR JACOBI CALCULATION JAC_EPS = 1.D-8 !** SET THE INITIAL GUESS x=0 !** SET INITIAL VALUES DO I = 1, M FVEC (I) = 0.D0 DO J = 1, N FJAC (I, J) = 0.D0 END DO END DO !** INITIALIZE SOLVER (ALLOCATE MEMORY, SET INITIAL VALUES) !** HANDLE IN/OUT: TR SOLVER HANDLE !** N IN: NUMBER OF FUNCTION VARIABLES !** M IN: DIMENSION OF FUNCTION VALUE !** X IN: SOLUTION VECTOR. CONTAINS VALUES X FOR F(X) !** EPS IN: PRECISIONS FOR STOP-CRITERIA !** ITER1 IN: MAXIMUM NUMBER OF ITERATIONS !** ITER2 IN: MAXIMUM NUMBER OF ITERATIONS OF CALCULATION OF TRIAL-STEP !** RS IN: INITIAL STEP BOUND IF (DTRNLSP_INIT (HANDLE, N, M, X, EPS, ITER1, ITER2, RS) & .NE. TR_SUCCESS) THEN !** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE PRINT *, '| ERROR IN DTRNLSP_INIT' !** RELEASE INTERNAL Intel(R) MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS. !** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS !** UNLESS YOU DISABLE Intel(R) MKL MEMORY MANAGER CALL MKL_FREE_BUFFERS !** AND STOP STOP 1 END IF !** CHECKS THE CORRECTNESS OF HANDLE AND ARRAYS CONTAINING JACOBIAN MATRIX, !** OBJECTIVE FUNCTION, LOWER AND UPPER BOUNDS, AND STOPPING CRITERIA. IF (DTRNLSP_CHECK (HANDLE, N, M, FJAC, FVEC, EPS, & INFO) .NE. TR_SUCCESS) THEN !** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE PRINT *, '| ERROR IN DTRNLSPBC_INIT' !** RELEASE INTERNAL Intel(R) MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS. !** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS !** UNLESS YOU DISABLE Intel(R) MKL MEMORY MANAGER CALL MKL_FREE_BUFFERS !** AND STOP STOP 1 ELSE !** THE HANDLE IS NOT VALID. IF( INFO(1) .NE. 0 .OR. !** THE FJAC ARRAY IS NOT VALID. & INFO(2) .NE. 0 .OR. !** THE FVEC ARRAY IS NOT VALID. & INFO(3) .NE. 0 .OR. !** THE EPS ARRAY IS NOT VALID. & INFO(4) .NE. 0 ) THEN PRINT *, '| INPUT PARAMETERS ARE NOT VALID' !** RELEASE INTERNAL Intel(R) MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS. !** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS !** UNLESS YOU DISABLE Intel(R) MKL MEMORY MANAGER CALL MKL_FREE_BUFFERS !** AND STOP STOP 1 END IF END IF !** SET INITIAL RCI CYCLE VARIABLES RCI_REQUEST = 0 SUCCESSFUL = 0 !** RCI CYCLE DO WHILE (SUCCESSFUL == 0) !** CALL TR SOLVER !** HANDLE IN/OUT: TR SOLVER HANDLE !** FVEC IN: VECTOR !** FJAC IN: JACOBI MATRIX !** RCI_REQUEST IN/OUT: RETURN NUMBER WHICH DENOTE NEXT STEP FOR PERFORMING IF (DTRNLSP_SOLVE (HANDLE, FVEC, FJAC, RCI_REQUEST) & .NE. TR_SUCCESS) THEN !** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE PRINT *, '| ERROR IN DTRNLSP_SOLVE' !** RELEASE INTERNAL Intel(R) MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS. !** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS !** UNLESS YOU DISABLE Intel(R) MKL MEMORY MANAGER CALL MKL_FREE_BUFFERS !** AND STOP STOP 1 END IF !** ACCORDING WITH RCI_REQUEST VALUE WE DO NEXT STEP SELECT CASE (RCI_REQUEST) CASE (-1, -2, -3, -4, -5, -6) !** STOP RCI CYCLE SUCCESSFUL = 1 CASE (1) !** RECALCULATE FUNCTION VALUE !** M IN: DIMENSION OF FUNCTION VALUE !** N IN: NUMBER OF FUNCTION VARIABLES !** X IN: SOLUTION VECTOR !** FVEC OUT: FUNCTION VALUE F(X) k=190000/3/(1-2*0.3); g=190000/(1+0.3)/2; pe= -167.6485 qe= 500.0385; CALL EXTENDED_POWELL (M, N, X, FVEC) CASE (2) !** COMPUTE JACOBI MATRIX !** EXTENDED_POWELL IN: EXTERNAL OBJECTIVE FUNCTION !** N IN: NUMBER OF FUNCTION VARIABLES !** M IN: DIMENSION OF FUNCTION VALUE !** FJAC OUT: JACOBI MATRIX !** X IN: SOLUTION VECTOR !** JAC_EPS IN: JACOBI CALCULATION PRECISION IF (DJACOBI (EXTENDED_POWELL, N, M, FJAC, X, JAC_EPS) & .NE. TR_SUCCESS) THEN !** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE PRINT *, '| ERROR IN DJACOBI' !** RELEASE INTERNAL Intel(R) MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS. !** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS !** UNLESS YOU DISABLE Intel(R) MKL MEMORY MANAGER CALL MKL_FREE_BUFFERS !** AND STOP STOP 1 END IF ENDSELECT END DO !** GET SOLUTION STATUSES !** HANDLE IN: TR SOLVER HANDLE !** ITER OUT: NUMBER OF ITERATIONS !** ST_CR OUT: NUMBER OF STOP CRITERION !** R1 OUT: INITIAL RESIDUALS !** R2 OUT: FINAL RESIDUALS IF (DTRNLSP_GET (HANDLE, ITER, ST_CR, R1, R2) & .NE. TR_SUCCESS) THEN !** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE PRINT *, '| ERROR IN DTRNLSP_GET' !** RELEASE INTERNAL Intel(R) MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS. !** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS !** UNLESS YOU DISABLE Intel(R) MKL MEMORY MANAGER CALL MKL_FREE_BUFFERS !** AND STOP STOP 1 END IF !** FREE HANDLE MEMORY IF (DTRNLSP_DELETE (HANDLE) .NE. TR_SUCCESS) THEN !** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE PRINT *, '| ERROR IN DTRNLSP_DELETE' !** RELEASE INTERNAL Intel(R) MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS. !** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS !** UNLESS YOU DISABLE Intel(R) MKL MEMORY MANAGER CALL MKL_FREE_BUFFERS !** AND STOP STOP 1 END IF !** RELEASE INTERNAL Intel(R) MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS. !** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS !** UNLESS YOU DISABLE Intel(R) MKL MEMORY MANAGER CALL MKL_FREE_BUFFERS !** IF FINAL RESIDUAL LESS THEN REQUIRED PRECISION THEN PRINT PASS IF (R2 .LT. 1.D-6) THEN PRINT *, '| DTRNLSP POWELL............PASS' print *,"x=",x print *,"r2=",r2 !pause !STOP 0 !** ELSE PRINT FAILED ELSE PRINT *, '| DTRNLSP POWELL............FAILED' print *, 'r2=',r2 print *, 'x=',x !STOP 1 END IF END PROGRAM EXAMPLE_DTRNLSP_POWELL !** ROUTINE FOR EXTENDED POWELL FUNCTION CALCULATION !** M IN: DIMENSION OF FUNCTION VALUE !** N IN: NUMBER OF FUNCTION VARIABLES !** X IN: VECTOR FOR FUNCTION CALCULATING !** F OUT: FUNCTION VALUE F(X) SUBROUTINE EXTENDED_POWELL (M, N, X, F) IMPLICIT NONE INTEGER M, N DOUBLE PRECISION X(*),F(*),K,G,PE,QE,SM & ,P,Q,PEEQ,FREAL,DS,DC,A INTEGER I COMMON PE,QE,K,G P=PE+K*X(1); Q=QE-3*G*X(2); PEEQ=X(3); SM=500+1000*(PEEQ); FREAL=X(4); DS=SINH(-1.5*P/SM); DC=COSH(-1.5*P/SM); A=0.04*EXP(-0.5*((PEEQ-0.3)/0.1)**2)/((2*3.14159)**0.5)/0.1; F(1) = X(1)* 2*Q/SM/SM+X(2)*(-3*1.5*FREAL*DS/SM); F(2) = (Q/SM)**2+2*FREAL*1.5*DC-(1+2.25*FREAL**2); F(3) = X(3)-(-P*X(1)+Q*X(2))/(1-FREAL)/SM F(4) = X(4)-(1-FREAL)*X(1)-A*PEEQ END SUBROUTINE EXTENDED_POWELL
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Dear customer,
The calculation relay on EPS, you set the EPS as 1E-5, it is too large for your sample. You'd better set smaller than the smallest value of x you would like to get, for instance 1.D-17. The result of MKL should be same as Matlab. That means all result of x would not be smaller than 1.D-17, thus the function value will also be different. You could print the result of function "FVEC" to check with it. Thanks.
Best regards,
Fiona
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Dear customer,
In MKL the final result is 1.540586674360078E-004, failed because you set the limitation for PASS is 1.D-6. I tried with your matlab code, seems there's an error. Could you please provide the output for matlab program and the init value & result value of x(1)~x(4) and the final f(x) value. Thanks.
Best regards,
Fiona
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
A good rule to follow when solving nonlinear equations is not to use all zeros as the guess for the solution, and a second rule is not to give up when an arbitrary initial guess fails to yield a solution.
I used x0 = [1d-7, 1d-7, 1d-7, 1d-7], and the solution found by TRNLSP was x = [3.3209d-12, 1.7481d-07, 1.7482d-07, 3.1324d-10] .
Perhaps you need to re-scale your variables x and F so that you can work with reasonable values, or use only a relative error in the termination criterion.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Dear Fiona, Thank you for your prompt reply. I just changed "x0=0" to "x0=[0,0,0,0]" in the Matlab code; Here is a table showing the input (x0)j (j=1,2,3,4), the solution x, and the residual F. The correct code is posted below the table. Seems like Matlab can get the maxinum residual to a order of 1e-14.
index j | 1 | 2 | 3 | 4 |
initial guess x0 | 0 | 0 | 0 | 0 |
Final value of x | 6.392630365754150e-17 | 1.74817817657479e-07 | 1.74817817600321e-07 | 3.09908247048591e-10 |
Final value of F | -4.26428382337290e-28 | -1.78745906964650e-14 | -4.10047855672419e-17 | -4.78885875879918e-21 |
!*********save the function as EQs.m file*****
function F = EQs(x)
k=190000/3/(1-2*0.3); g=190000/(1+0.3)/2; pe=-167.6485 ; qe=500.0385 ; p=pe+k*x(1); q=qe-3*g*x(2); peeq=x(3); sm=500+1000*(peeq); freal=x(4); ds=sinh(-1.5*p/sm); dc=cosh(-1.5*p/sm); a=0.04*exp(-0.5*((peeq-0.3)/0.1)^2)/((2*3.14159)^0.5)/0.1; F(1) = x(1)* 2*q/sm/sm+x(2)*(-3*1.5*freal*ds/sm); F(2) = (q/sm)^2+2*freal*1.5*dc-(1+2.25*freal^2); F(3)=x(3)-(-p*x(1)+q*x(2))/(1-freal)/sm F(4)=x(4)-(1-freal)*x(1)-a*peeq
!**************end of EQs.m file**************
>>options = optimoptions('fsolve','Display','iter','PlotFcn',@optimplotfirstorderopt); >>x0=[0,0,0,0]; >>[x,fval] = fsolve(@EQs,x0,options)
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Hi mecej4,
Thanks for your reply, your result is very close to what I have got via Matlab using initial guess x0=[0,0,0,0], pls see the result in my above comment.
- However, I tried x0=[1d-7,1d-7,1d-7,1d-7,], failed still. Have you changed other line of Fortran code?
- What do you mean by "rescale" x and F? Could you pls elaborate on how to do so in TRNLSP ?
Kind regards,
P
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Dear customer,
The calculation relay on EPS, you set the EPS as 1E-5, it is too large for your sample. You'd better set smaller than the smallest value of x you would like to get, for instance 1.D-17. The result of MKL should be same as Matlab. That means all result of x would not be smaller than 1.D-17, thus the function value will also be different. You could print the result of function "FVEC" to check with it. Thanks.
Best regards,
Fiona

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page