CODEWITHSUNDEEP
matlabnewtons-method
Miro Feliciano
Miro Feliciano

Reputation: 79

Is there a faster way to solve my newton algorithm?

I have this Newton algorithm which seems to be taking forever to find a solution. To the point where i let it sit for a whole night and nothing happend. I was wondering if it is something I put wrong or maybe something I am missing.

Gobs=[0.61 1.14 2.33 4.76 6.65 4.77 2.38 1.13 0.59];
x=[-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.5];
syms m h
p1=[m;h];

F = ((((6.67*p1(1)*p1(2))/((x(1)^2)+p1(2)^2)^(3/2))-Gobs(1))^2+(((6.67*p1(1)*p1(2))/((x(2)^2)+p1(2)^2)^(3/2))-Gobs(2))^2 + (((6.67*p1(1)*p1(2))/((x(3)^2)+p1(2)^2)^(3/2))-Gobs(3))^2 + (((6.67*p1(1)*p1(2))/((x(4)^2)+p1(2)^2)^(3/2))-Gobs(4))^2 + (((6.67*p1(1)*p1(2))/((x(5)^2)+p1(2)^2)^(3/2))-Gobs(5))^2 + (((6.67*p1(1)*p1(2))/((x(6)^2)+p1(2)^2)^(3/2))-Gobs(6))^2 + (((6.67*p1(1)*p1(2))/((x(7)^2)+p1(2)^2)^(3/2))-Gobs(7))^2 + (((6.67*p1(1)*p1(2))/((x(8)^2)+p1(2)^2)^(3/2))-Gobs(8))^2 + (((6.67*p1(1)*p1(2))/((x(9)^2)+p1(2)^2)^(3/2))-Gobs(9))^2);

dfm = diff(F,'m');
dfh = diff(F,'h');
d2fm = diff(dfm,'m'); 
d2fh = diff(dfh,'h');
dfmh = diff(dfm,'h');
dfhm = diff(dfh,'m');
G = [dfm;dfh];
Hnewton = [d2fm dfmh; dfmh d2fh];    

m =1.6;
h =1.7;

p1 = subs(p1);
n=1;
tm(n)=p1(1);
th(n)=p1(2);
for i=1:2
F1 = subs(F);
p2 = p1-(subs(Hnewton)+ (100/n)*eye(2))\subs(G);
m=p2(1);
p=p2(2);
F2 = subs(F);
if (F2)>0.9;
p1 = p2;
m=p2(1);
h=p2(2);
else break
end
n=n+1;
tm(n)=p2(1);
th(n)=p2(2);
end

plot(th,tm,'r-','LineWidth',1.4)

Upvotes: 1

Views: 126

Answers (1)

Adam
Adam

Reputation: 2777

Define two versions of the function to be optimized:

  • F_eval function handle type used for evaluating F,

  • F syms function used for the gradient and the hessian matrix

To evaluate the gradient G and the hessian matrix Hnewton, used their corresponding function handle

% Gradient
G_eval = matlabFunction(G,'Vars',{[m h]});

% Hessian 
Hnewton_eval = matlabFunction(Hnewtonn,'Vars',{[m h]});

Also the loop stopping condition is wrong:

use the difference between two successive function evaluation

Read through the code below 


clc
clear
syms m h
p1 = [m h];

Gobs=[0.61 1.14 2.33 4.76 6.65 4.77 2.38 1.13 0.59];
x=[-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.5];
% Define the input known array x first 
% Used only for gradient and hessian, syms function  
F =(...
    (((6.67*p1(1)*p1(2))/((x(1)^2)+p1(2)^2)^(3/2))-Gobs(1))^2 ...
    +(((6.67*p1(1)*p1(2))/((x(2)^2)+p1(2)^2)^(3/2))-Gobs(2))^2 ...
    +(((6.67*p1(1)*p1(2))/((x(3)^2)+p1(2)^2)^(3/2))-Gobs(3))^2 ...
    +(((6.67*p1(1)*p1(2))/((x(4)^2)+p1(2)^2)^(3/2))-Gobs(4))^2 ... 
    +(((6.67*p1(1)*p1(2))/((x(5)^2)+p1(2)^2)^(3/2))-Gobs(5))^2 ... 
    +(((6.67*p1(1)*p1(2))/((x(6)^2)+p1(2)^2)^(3/2))-Gobs(6))^2 ...
    +(((6.67*p1(1)*p1(2))/((x(7)^2)+p1(2)^2)^(3/2))-Gobs(7))^2 ...
    +(((6.67*p1(1)*p1(2))/((x(8)^2)+p1(2)^2)^(3/2))-Gobs(8))^2 ... 
    +(((6.67*p1(1)*p1(2))/((x(9)^2)+p1(2)^2)^(3/2))-Gobs(9))^2 ...
   );
% Used for evaluating the function , function handle 
F_eval = @(p)...
    (...
    (((6.67*p(1)*p(2))/((x(1)^2)+p(2)^2)^(3/2))-Gobs(1))^2 ...
    +(((6.67*p(1)*p(2))/((x(2)^2)+p(2)^2)^(3/2))-Gobs(2))^2 ...
    +(((6.67*p(1)*p(2))/((x(3)^2)+p(2)^2)^(3/2))-Gobs(3))^2 ...
    +(((6.67*p(1)*p(2))/((x(4)^2)+p(2)^2)^(3/2))-Gobs(4))^2 ... 
    +(((6.67*p(1)*p(2))/((x(5)^2)+p(2)^2)^(3/2))-Gobs(5))^2 ... 
    +(((6.67*p(1)*p(2))/((x(6)^2)+p(2)^2)^(3/2))-Gobs(6))^2 ...
    +(((6.67*p(1)*p(2))/((x(7)^2)+p(2)^2)^(3/2))-Gobs(7))^2 ...
    +(((6.67*p(1)*p(2))/((x(8)^2)+p(2)^2)^(3/2))-Gobs(8))^2 ... 
    +(((6.67*p(1)*p(2))/((x(9)^2)+p(2)^2)^(3/2))-Gobs(9))^2 ...
   );


% Gradient and Hessian in syms 
G = gradient(F);
Hnewton = hessian(F);

%% Tranform into function handle


G_eval = matlabFunction(G,'Vars',{[m h]});

Hnewton_eval = matlabFunction(Hnewton,'Vars',{[m h]});

%% Starting guess
m = 1.6;
h =1.7;
n = 1;
pold = [m h];

%% Record 
% First column is m, second is h
solution_history = pold;

%% Infinite loop
while true

     % Use function handle only for evaluation
     F1 = F_eval(pold);

     psol = pold -(Hnewton_eval(pold)+ (100/n)*eye(2))\G_eval(pold);

     [l,~] = size(psol);

     F2_best = F_eval(psol(1, :));
     p_best = psol(1, :);

     % Find the solution which minimize the most F
     for i = 1:l
         if  F_eval(psol(i, :)) < F2_best
             F2_best = F_eval(psol(i, :));
             p_best = psol(1, :);
         end
     end

    F2 = F2_best;
    pnew = p_best;
     solution_history = [solution_history; pnew];

     n = n + 1;

     pold = pnew;
     % Insert stopping condition here  
      if n == 5
          break;
      end

end

tm = solution_history(:, 1);
th = solution_history(:, 2);

plot(th,tm,'r-','LineWidth',1.4)

Upvotes: 1

Related Questions