1

I am trying to find a solution for a non-linear system in MATLAB using the fsolve function. I have an idea of the region the solution is, so my initial points are a random variation inside this area. The code follows:

%MATLAB parameters
digits = 32;
format shortEng
%Initialization of solving process
min_norm = realmax;
options = optimoptions('fsolve', 'Algorithm', 'trust-region-dogleg',... 
    'Display', 'off', 'FunValCheck', 'on', 'Jacobian', 'on',...
    'DerivativeCheck', 'on');
%Solving loop
while 1
    %Refreshes seed after 100000 iterations
    rng('shuffle');
    for n=1:100000
        %Initial points are randomly placed in solving region
        x_0 = zeros(1, 3);
        x_0(1) = 20*(10^-3)+ abs(3*(10^-3)*randn);
        x_0(2) = abs(10^(-90) + 10^-89*randn);
        x_0(3) = abs(10*randn);
        %Solving
        [x, fval] = fsolve(@fbnd, x_0, options);
        norm_fval = norm(fval);
        %If norm is minimum, result is displayed
        if all(x > 0.0) && (norm_fval < min_norm)
            iter_solu = x;
            display(norm_fval);
            display(iter_solu);
            min_norm = norm_fval;
        end
    end
end

The function to be optimized and its jacobian follow:

function [F, J] = fbnd(x)
    F(1) = x(1) - x(2)*(exp(7.56/(1.5*0.025))-1);
    F(2) = x(1) - x(2)*(exp((x(3)*0.02)/(1.5*0.025))-1)-0.02;
    F(3) = x(1) - x(2)*(exp((6.06+x(3)*0.018)/(1.5*0.025))-1)-0.018;
    J = [[1.0, -3.5790482371355382991651020082214*10^87, 0.0];...
         [1.0, 1.0-exp((8/15)*x(3)),-(8/15)*x(2)*exp((8/15)*x(3))];...
         [1.0, 1.0-exp(0.48*x(3)+161.6), -0.48*x(2)*exp(0.48*x(3)+161.6)]];
end

When I turn DerivativeCheck on, I get the following error:

Objective function derivatives:
Maximum relative difference between user-supplied 
and finite-difference derivatives = 1.
  User-supplied derivative element (1,1):     1
  Finite-difference derivative element (1,1): 0

For some reason, it things dF(1)/dx(1) should be zero instead of 1, when it is obvious it is 1. So I went ahead and replaced the J(1,1) by 0, and now I get the same problem with J(3,1):

Objective function derivatives:
Maximum relative difference between user-supplied 
and finite-difference derivatives = 1.
  User-supplied derivative element (3,1):     1
  Finite-difference derivative element (3,1): 0

So, I replaced J(3,1) by 0 and it catches only minor precision errors. Is there an explanation why it thinks J(1,1) and J(3,1) should be zero? It seems really strange to me.

πάντα ῥεῖ
  • 1
  • 13
  • 116
  • 190
gstorto
  • 165
  • 8
  • 3
    It is only a supposition, but since your Jacobian have extremely different scales the finite differences might not be accurate enough. Moreover, `validateFirstDerivatives` has an internal component-wise relative tolerance of 1e-6. I suggest you to scale your problem, so that F_1 and F_3 have the same magnitude of F_2. – Jommy Jul 24 '15 at 09:17
  • 2
    Well with elements ranging from `1.0` to `-3.58*10^87` it's little surprise you are encountering problems. `rcond(J)` is `2.7940e-88`. The function itself is numerically badly scaled as well. If you compute `eps((exp(7.56/(1.5*0.025))-1))` you get `4.4171e+71` which says that finite forward differencing is likely to give inaccurate or incorrect results. – transversality condition Jul 24 '15 at 09:45
  • @Jommy thanks for the comments. Yeah. After I posted I realized it might have some connection to the different magnitudes of the problem. I will read the article Matlab has about scaling. I think it will cause problems not only on the Jacobian but also in the solving as well. – gstorto Jul 24 '15 at 13:33

0 Answers0