MATLAB: What did starting points could be with the MATLAB script

starting values

The following MATLAB script;

    function [x,fx,xx] = newtons(f,x0,TolX,MaxIter,varargin)
    %newtons.m to solve a set of nonlinear eqs f1(x)=0, f2(x)=0, f3(x)=0,
    %f4(x)=0,..
    %input:  f   = 1^st-order vector ftn equivalent to a set of equations
    %        x0  = the initial guess of the solution
    %       TolX = the upper limit of |x(k) - x(k-1)|
    %    MaxIter = the maximum # of iteration
    %output: x   = the point which the algorithm has reached
    %        fx  = f(x(last))
    %        xx  = the history of x
    h = 1e-4; TolFun = eps; EPS = 1e-6;
    fx = feval(f,x0,varargin{:});
    Nf = length(fx); Nx = length(x0);
    if Nf ~= Nx, error('Incompatible dimensions of f and x0!'); end
    if nargin < 4, MaxIter = 200; end
    if nargin < 3, TolX = EPS; end
    xx(1,:) = x0(:).'; %Initialize the solution as the initial row vector
    %fx0 = norm(fx);                              %(1)
    for k = 1: MaxIter
      dx = -jacob(f,xx(k,:),h,varargin{:})\fx(:);%-[dfdx]^-1*fx
      %for I = 1; 3 %damping to avoid divergence %(2)
      %dx = dx/2;                                %(3)
      xx(k + 1,:) = xx(k,:) + dx.';
      fx = feval(f,xx(k + 1,:),varargin{:}); fxn = norm(fx);
      % if fxn < fx0, break; end                 %(4)
      %end                                       %(5)
      if fxn < TolFun | norm(dx) < TolX, break; end
      %fx0 = fxn;                                %(6)
    end
    x = xx(k + 1,:);
    if k == MaxIter, fprintf('The best in %d iterations\n',MaxIter), end
    function g = jacob(f,x,h,varargin) %Jacobian of f(x)
    if nargin < 3, h = 1e-6; end
    h2 = 2*h; N = length(x); x = x(:).'; I = eye(N);
    for n = 1:N
      g(:,n) = (feval(f,x + I(n,:)*h,varargin{:}) ...
                   -feval(f,x - I(n,:)*h,varargin{:}))'/h2;
    end
was to solve the following system of nonlinear equations given as;
    y(1) = 5-x(1)*(6^x(2)*exp(6*x(3))+12^x(2)*exp(12*x(3))+18^x(2)*exp(18*x(3))+24^x(2)*exp(24*x(3))+30^x(2)*exp(30*x(3)))-2*x(1)*x(4)*(6^x(2)*exp(-x(1)*6^x(2)*exp(6*x(3))+(6*x(3)))*(1-x(4)*exp(-x(1)*6^x(2)*exp(6*x(3))))^(-1)+12^x(2)*exp(-x(1)*12^x(2)*exp(12*x(3))+(12*x(3)))*(1-x(4)*exp(-x(1)*12^x(2)*exp(12*x(3))))^(-1)+18^x(2)*exp(-x(1)*18^x(2)*exp(18*x(3))+(18*x(3)))*(1-x(4)*exp(-x(1)*18^x(2)*exp(18*x(3))))^(-1)+24^x(2)*exp(-x(1)*24^x(2)*exp(24*x(3))+(24*x(3)))*(1-x(4)*exp(-x(1)*24^x(2)*exp(24*x(3))))^(-1)+30^x(2)*exp(-x(1)*30^x(2)*exp(30*x(3))+(30*x(3)))*(1-x(4)*exp(-x(1)*30^x(2)*exp(30*x(3))))^(-1))=0;
    y(2) = (log(6)+log(12)+log(18)+log(24)+log(30))+((x(2)+6*x(3))^(-1)+(x(2)+12*x(3))^(-1)+(x(2)+18*x(3))^(-1)+(x(2)+24*x(3))^(-1)+(x(2)+30*x(3))^(-1))-x(1)*(6^x(2)*exp(6*x(3))*log(6)+12^x(2)*exp(12*x(3))*log(12)+18^x(2)*exp(18*x(3))*log(18)+24^x(2)*exp(24*x(3))*log(24)+30^x(2)*exp(30*x(3))*log(30))-2*x(2)*x(4)*(6^x(2)*exp(-x(1)*6^x(2)*exp(6*x(3))+(6*x(3)))*log(6)*(1-x(4)*exp(-x(1)*6^x(2)*exp(6*x(3))))^(-1)+12^x(2)*exp(-x(1)*12^x(2)*exp(12*x(3))+(12*x(3)))*log(12)*(1-x(4)*exp(-x(1)*12^x(2)*exp(12*x(3))))^(-1)+18^x(2)*exp(-x(1)*18^x(2)*exp(18*x(3))+(18*x(3)))*log(18)*(1-x(4)*exp(-x(1)*18^x(2)*exp(18*x(3))))^(-1)+24^x(2)*exp(-x(1)*24^x(2)*exp(24*x(3))+(24*x(3)))*log(24)*(1-x(4)*exp(-x(1)*24^x(2)*exp(24*x(3))))^(-1)+30^x(2)*exp(-x(1)*30^x(2)*exp(30*x(3))+(30*x(3)))*log(30)*(1-x(4)*exp(-x(1)*30^x(2)*exp(30*x(3))))^(-1))=0;
    y(3) = (6*(x(2)+6*x(3))^(-1)+12*(x(2)+12*x(3))^(-1)+18*(x(2)+18*x(3))^(-1)+24*(x(2)+24*x(3))^(-1)+30*(x(2)+30*x(3))^(-1))-x(1)*(6^(x(2)+1)*exp(6*x(3))+12^(x(2)+1)*exp(12*x(3))+18^(x(2)+1)*exp(18*x(3))+24^(x(2)+1)*exp(24*x(3))+30^(x(2)+1)*exp(30*x(3)))-(6+12+18+24+30)-2*x(1)*x(4)*(6^(x(2)+1)*exp(-x(1)*6^x(2)*exp(6*x(3))+(6*x(3)))*(1-x(4)*exp(-x(1)*6^x(2)*exp(6*x(3))))^(-1)+12^(x(2)+1)*exp(-x(1)*12^x(2)*exp(12*x(3))+(12*x(3)))*(1-x(4)*exp(-x(1)*12^x(2)*exp(12*x(3))))^(-1)+18^(x(2)+1)*exp(-x(1)*18^x(2)*exp(18*x(3))+(18*x(3)))*(1-x(4)*exp(-x(1)*18^x(2)*exp(18*x(3))))^(-1)+24^(x(2)+1)*exp(-x(1)*24^x(2)*exp(24*x(3))+(24*x(3)))*(1-x(4)*exp(-x(1)*24^x(2)*exp(24*x(3))))^(-1)+30^(x(2)+1)*exp(-x(1)*30^x(2)*exp(30*x(3))+(30*x(3)))*(1-x(4)*exp(-x(1)*30^x(2)*exp(30*x(3))))^(-1))=0;
    y(4) = 5*(1-x(4))^(-1)-2*(exp(-x(1)*6^x(2)*exp(6*x(3)))*(1-x(4)*exp(-x(1)*6^x(2)*exp(6*x(3))))^(-1)+exp(-x(1)*12^x(2)*exp(12*x(3)))*(1-x(4)*exp(-x(1)*12^x(2)*exp(12*x(3))))^(-1)+exp(-x(1)*18^x(2)*exp(18*x(3)))*(1-x(4)*exp(-x(1)*18^x(2)*exp(18*x(3))))^(-1)+exp(-x(1)*24^x(2)*exp(24*x(3)))*(1-x(4)*exp(-x(1)*24^x(2)*exp(24*x(3))))^(-1)+exp(-x(1)*30^x(2)*exp(30*x(3)))*(1-x(4)*exp(-x(1)*30^x(2)*exp(30*x(3))))^(-1))=0;
    end

I have tried various starting values for x0 but the output is always NaN NaN NaN NaN.Please, i need help for the starting values. Thanks.

Best Answer

You want us to give you starting values for some arbitrary set of highly nonlinear equations, to make Newton's method converge? Sorry, I should have said to divine them? The crystal ball is so foggy. Sorry. Not gonna happen.

Newton's method is unconstrained. So it can trivially diverge into outer space, especially with such a nasty set of nonlinear equations.

So, first, learn to use a better solver instead of homegrown code. One that allows constraints.

Second, try to come up with a set of intelligent constraints.

Third, spend some time investigating the parameter space. At your starting values, are these equations returning NaNs, probably because your equations are so nasty? Exponents and exponentials are bad things. If you don't appreciate that, then it is time to do some reading about floating point arithmetic. Learn about the dynamic range of a double precision number.

Ok, maybe you don't believe me about dynamic range, or don't appreciate what it means?

Look at this code fragment from just ONE of your equations:

exp(-x(1)*18^x(2)*exp(18*x(3)))

So you try to compute exp of the expression:

-x(1)*18^x(2)*exp(18*x(3))

Just what values of that product do you think will not cause exp to overflow or underflow? For example, depending on the value of x(3), exp(18*x(3)) will either be immense, overflowing to inf, or underflow to zero.

How big can x(3) be? Not very big.

exp(709)
ans =
     8.2184e+307
exp(710)
ans =
     Inf

So, if x(3) is larger than roughly 39, exp(18*x(3)) will overflow.

709/18
ans =
       39.389

exp will as quickly return zero if x(3) is a negative number too. Don't forget that Newton's method has no constraints.

Then, you multiply the result of exp(18*x(3)) with terms like 18^x(2).

How large does x(2) need to be in order for that term to overflow? Again, not that large. Then, after you multiply those two results by -x(1), you try to compute exp of THAT mess?

So again, all of these terms combine to make something that is simply insolvable in floating point arithmetic, when done in a reasonable number of digits of precision. It is especially so when you have NO idea what values of those parameters will yield finite values for those equations.

Worse, given the utterly incredibly highly nonlinear equations that you have posed, you honestly intend to use a tool that computes a Jacobian using simple finite differences with a fixed step size? Sorry, but it is simply incredible to hope these computations will yield anything meaningful.

Yes, I know, you won't believe me. I'll define those functions, then see if I can even evaluate them in symbolic form.

y{1} = @(x) 5-x(1)*(6^x(2)*exp(6*x(3))+12^x(2)*exp(12*x(3))+18^x(2)*exp(18*x(3))+24^x(2)*exp(24*x(3))+30^x(2)*exp(30*x(3)))-2*x(1)*x(4)*(6^x(2)*exp(-x(1)*6^x(2)*exp(6*x(3))+(6*x(3)))*(1-x(4)*exp(-x(1)*6^x(2)*exp(6*x(3))))^(-1)+12^x(2)*exp(-x(1)*12^x(2)*exp(12*x(3))+(12*x(3)))*(1-x(4)*exp(-x(1)*12^x(2)*exp(12*x(3))))^(-1)+18^x(2)*exp(-x(1)*18^x(2)*exp(18*x(3))+(18*x(3)))*(1-x(4)*exp(-x(1)*18^x(2)*exp(18*x(3))))^(-1)+24^x(2)*exp(-x(1)*24^x(2)*exp(24*x(3))+(24*x(3)))*(1-x(4)*exp(-x(1)*24^x(2)*exp(24*x(3))))^(-1)+30^x(2)*exp(-x(1)*30^x(2)*exp(30*x(3))+(30*x(3)))*(1-x(4)*exp(-x(1)*30^x(2)*exp(30*x(3))))^(-1));
y{2} = @(x)  (log(6)+log(12)+log(18)+log(24)+log(30))+((x(2)+6*x(3))^(-1)+(x(2)+12*x(3))^(-1)+(x(2)+18*x(3))^(-1)+(x(2)+24*x(3))^(-1)+(x(2)+30*x(3))^(-1))-x(1)*(6^x(2)*exp(6*x(3))*log(6)+12^x(2)*exp(12*x(3))*log(12)+18^x(2)*exp(18*x(3))*log(18)+24^x(2)*exp(24*x(3))*log(24)+30^x(2)*exp(30*x(3))*log(30))-2*x(2)*x(4)*(6^x(2)*exp(-x(1)*6^x(2)*exp(6*x(3))+(6*x(3)))*log(6)*(1-x(4)*exp(-x(1)*6^x(2)*exp(6*x(3))))^(-1)+12^x(2)*exp(-x(1)*12^x(2)*exp(12*x(3))+(12*x(3)))*log(12)*(1-x(4)*exp(-x(1)*12^x(2)*exp(12*x(3))))^(-1)+18^x(2)*exp(-x(1)*18^x(2)*exp(18*x(3))+(18*x(3)))*log(18)*(1-x(4)*exp(-x(1)*18^x(2)*exp(18*x(3))))^(-1)+24^x(2)*exp(-x(1)*24^x(2)*exp(24*x(3))+(24*x(3)))*log(24)*(1-x(4)*exp(-x(1)*24^x(2)*exp(24*x(3))))^(-1)+30^x(2)*exp(-x(1)*30^x(2)*exp(30*x(3))+(30*x(3)))*log(30)*(1-x(4)*exp(-x(1)*30^x(2)*exp(30*x(3))))^(-1));
y{3} = @(x) (6*(x(2)+6*x(3))^(-1)+12*(x(2)+12*x(3))^(-1)+18*(x(2)+18*x(3))^(-1)+24*(x(2)+24*x(3))^(-1)+30*(x(2)+30*x(3))^(-1))-x(1)*(6^(x(2)+1)*exp(6*x(3))+12^(x(2)+1)*exp(12*x(3))+18^(x(2)+1)*exp(18*x(3))+24^(x(2)+1)*exp(24*x(3))+30^(x(2)+1)*exp(30*x(3)))-(6+12+18+24+30)-2*x(1)*x(4)*(6^(x(2)+1)*exp(-x(1)*6^x(2)*exp(6*x(3))+(6*x(3)))*(1-x(4)*exp(-x(1)*6^x(2)*exp(6*x(3))))^(-1)+12^(x(2)+1)*exp(-x(1)*12^x(2)*exp(12*x(3))+(12*x(3)))*(1-x(4)*exp(-x(1)*12^x(2)*exp(12*x(3))))^(-1)+18^(x(2)+1)*exp(-x(1)*18^x(2)*exp(18*x(3))+(18*x(3)))*(1-x(4)*exp(-x(1)*18^x(2)*exp(18*x(3))))^(-1)+24^(x(2)+1)*exp(-x(1)*24^x(2)*exp(24*x(3))+(24*x(3)))*(1-x(4)*exp(-x(1)*24^x(2)*exp(24*x(3))))^(-1)+30^(x(2)+1)*exp(-x(1)*30^x(2)*exp(30*x(3))+(30*x(3)))*(1-x(4)*exp(-x(1)*30^x(2)*exp(30*x(3))))^(-1));
y{4} = @(x) 5*(1-x(4))^(-1)-2*(exp(-x(1)*6^x(2)*exp(6*x(3)))*(1-x(4)*exp(-x(1)*6^x(2)*exp(6*x(3))))^(-1)+exp(-x(1)*12^x(2)*exp(12*x(3)))*(1-x(4)*exp(-x(1)*12^x(2)*exp(12*x(3))))^(-1)+exp(-x(1)*18^x(2)*exp(18*x(3)))*(1-x(4)*exp(-x(1)*18^x(2)*exp(18*x(3))))^(-1)+exp(-x(1)*24^x(2)*exp(24*x(3)))*(1-x(4)*exp(-x(1)*24^x(2)*exp(24*x(3))))^(-1)+exp(-x(1)*30^x(2)*exp(30*x(3)))*(1-x(4)*exp(-x(1)*30^x(2)*exp(30*x(3))))^(-1));

I'll pick some arbitrary, reasonably small values for x.

x0 = sym([1 1 1 1]);
vpa(y{1}(x0))
ans =
-321231160211767.64948297735800268
vpa(y{2}(x0))
ans =
-1092428114449996.7550800164309852
vpa(y{3}(x0))
ans =
-9633106155000105.0129906212143782
vpa(y{4}(x0))
Error using symengine
Division by zero.
Error in sym/privBinaryOp (line 908)
            Csym = mupadmex(op,args{1}.s, args{2}.s, varargin{:});
Error in  ^  (line 266)
        B = privBinaryOp(A, p, 'symobj::mpower');
Error in
@(x)5*(1-x(4))^(-1)-2*(exp(-x(1)*6^x(2)*exp(6*x(3)))*(1-x(4)*exp(-x(1)*6^x(2)*exp(6*x(3))))^(-1)+exp(-x(1)*12^x(2)*exp(12*x(3)))*(1-x(4)*exp(-x(1)*12^x(2)*exp(12*x(3))))^(-1)+exp(-x(1)*18^x(2)*exp(18*x(3)))*(1-x(4)*exp(-x(1)*18^x(2)*exp(18*x(3))))^(-1)+exp(-x(1)*24^x(2)*exp(24*x(3)))*(1-x(4)*exp(-x(1)*24^x(2)*exp(24*x(3))))^(-1)+exp(-x(1)*30^x(2)*exp(30*x(3)))*(1-x(4)*exp(-x(1)*30^x(2)*exp(30*x(3))))^(-1))

As it turns out, smaller values of x, will generate smaller results. As I said, investigate the parameter space. You created this problem. You need to understand what is happening! However, when x is entirely zero, you will find more singularities, more divide by zeros. Newton's method will fail miserably here, unless you can find exceptionally good starting values. However, I did say to find a better optimizer than what appears to be code written by a student taking their first class in numerical methods.

So it looks like there is a divide by zero when x(4)==1. The point is, there are lots of singularities, various type of overflows and underflows in these equations, even if you use super-high precision. You won't get any useful answer from what you are trying to do though, certainly not as you are trying to solve the problem.

Computers are NOT infinitely powerful, they only seem so because we see them that way on tv or in the movies.

Best Answer

Related Solutions

MATLAB: Writing several sums of terms in matlab

MATLAB: Some of the values of the estimated parameters do not conform

Related Question