Newton method for solving a nonlinear system of equations!

Hello,

I am trying to write a function that solves a non linear system of equations. In my exercise I have 2,3 and 4 dimensional equation system. I am not allowed to use fsolve function. I also have a function called "Score4" which is suppose to calculate the score of my function ( should be around 1 if the code that I wrote is correct).

Below you will find the score and the code that I wrote.

Thank you in advance for your time.

function [score,error,details] = Score4()

%% Parameter

% 2 dimensional equation system
A{1} = @(x)[4*cos(x(1)),2*cos(x(1)+x(2));
    5*sin(x(1)),1*sin(x(1)+x(2))];
x0{1} = [0,0]';
y{1} = [2,3]';
y0{1} = [-1,1]';

% 3 dimensional equation system
A{2} = @(x)[1,0,0;0,cos(x(1)),-sin(x(2));0,sin(x(2))^2,cos(x(1))];
x0{2} = [0;0];
y0{2} = [1;2;3];
y{2} = [1;0;1];

% 4 dimensional transformation matrix with 5d input
A{3} = @(x)[cos(x(1)),sin(x(1)),0,-x(3);
-sin(x(1))*cos(x(4)),cos(x(1))*cos(x(4)),sin(x(4))*sin(x(5)),-x(2)*sin(x(4))*cos(x(5));
sin(x(4))*sin(x(1)),-cos(x(1))*sin(x(4)),cos(x(4))*sin(x(5)),-x(2)*cos(x(4))*cos(x(5));
0,0,0,1];
x0{3} = [0;0;0;0;0];
y0{3} = [1;3;2;1];
y{3} = [1;2;-1;1];


%% Score Berechnung
score = 0;
error = '';
details = struct();
try
    % Iteration ueber testcases
     for i = 1:3
        % Aufstellen des nullstellenproblems
        RootProblem = @(z) A{i}(z)*y0{i}-y{i};
        T = 0;
        % Mitteln ueber Zeit
		timer = tic();
        for j = 1:10
            x = nlinSolve(RootProblem,x0{i});   
        end
        T = toc(timer)/10;
        output = A{i}(x)*y0{i}-y{i};
        Err = norm(output);
        Sc = (Err+0.01)*(T+0.5);
        score = score+Sc;
        
        details.(['time',num2str(i)]) = T;
        details.(['error',num2str(i)]) = Err;
     end
     score = 0.017 / score;
catch err
    score = 0;
    error = err.message;
end

end

To solve this problem I was thinking to use Newton method but I am having a lot of problems. Can someone give me some suggestions how to improve my code and what am I doing wrong. My code is as follows:

function x = nlinSolve (g,x0)


 
select = @(g,r,c) g(r,c);
G1 = @(x)select(g(x), 1, 1);    % selecting the first row of g(x) 
G2 =  @(x)select(g(x), 2, 1);   % selecting the second row
G3 =  @(x)select(g(x), 3, 1);
G4 =  @(x)select(g(x), 4, 1);
G5 =  @(x)select(g(x), 5, 1);

x = x0;         % x0 initial guess which is given in the score function 
for m=1:10 

%calculating Jacobian

for k=1:numel(x)
    for n =1: numel(x)

        Jakobian(k,n)= Derivative(eval(['G',num2str(k)]),x,x(n),n);

    end
end


x = x-inv(Jakobian)*g(x);  % Newton formula to find roots

end
end

function J = Derivative(G,x,xm,m)

% the funciton that calculates the derivative using numerical method

e = 0.0001;

x1=x;
x2=x;

x1(m) = xm;

x2(m) = xm+e;

J = (G(x2)-G(x1))/e;


end

Now this code is clearly not working, I am messing up something. Also the score should be around 1 at the end.

If anyone can help me that would be more than appritiated.

[Rigon]

If someone can tell me how to calculate the Jacobian of this type of functions. For example

A{1} = @(x)[4*cos(x(1)),2*cos(x(1)+x(2));
    5*sin(x(1)),1*sin(x(1)+x(2))];

g = g=@(x) A{i}(x)*y0{i} - y{i};

How to calculate the Jacobian of the function g. That probably would be really helpful.

Thank you again!

Answer 1

[Matt J]

How to calculate the Jacobian of the function g. If Ja is the Jacobian of A{i}(:), then the Jacobian Jgof g will be

Jg=kron(y0{i}.',eye(2))*Ja;

[Rigon]

Thank you. How to calculate Ja in this case? With 4 separate functions in a matrix 2x2 inside a function handle, I don't think I ever have done it. Thank you again for your time!

[Matt J]

You should supply the analytical form for Ja in the same way you have supplied it for A

Ja=@(x)[  4*sin(x(1)),     0;
         -5*cos(x(1)),     0;
         -2*sin(x(1)+x(2)), -2*sin(x(1)+x(2));
          cos(x(1)+x(2)) ,  cos(x(1)+x(2))]

[Rigon]

When I try to run the code for Jg there is an error:

Operator '*' is not supported for operands of type 'function_handle'.

Probably I should try without function handle for Ja just as matrix.

[Matt J]

Yes. It should be,

Jg= kron(y0{i}.',eye(2))*Ja(x);

All quantities on both sides of the above are numeric.

[Rigon]

Thank you, but still when I calculate the Jg, I am not able to find the inverse of Jg which I need to find the roots. x = x-inv(Jg)*g(x); There is a warning: Matrix is singular to working precision!

[Matt J]

It is possible for the Jacobian to be singular, which is why pure Newton's method is seldom used. If y0{i}=[0;0] for example, the Jacobian will be Jg=zeros(2) for all x.

[Rigon]

Ok for example in my case I have:

A{1} = @(x)[4*cos(x(1)),2*cos(x(1)+x(2));
    5*sin(x(1)),1*sin(x(1)+x(2))];
x0{1} = [0,0]';
y{1} = [2,3]';
y0{1} = [-1,1]';
g= @(x) A{1}(x)*y0{1}-y{1};
x=x0{1};

for m=1:10

Ja=@(x) [  4*sin(x(1)),     0;
         -5*cos(x(1)),     0;
         -2*sin(x(1)+x(2)), -2*sin(x(1)+x(2));
          cos(x(1)+x(2)) ,  cos(x(1)+x(2))];

Jg=kron(y0{1}.',eye(2))*Ja(x);

x = x-inv(Jg)*g(x);

end

I can not find the x in this case and is supose to be the same as

x = fsolve(g,x)

which I can calculate very easy.

[Matt J]

A{1} = @(x)[4*cos(x(1)),2*cos(x(1)+x(2));
    5*sin(x(1)),1*sin(x(1)+x(2))];
x0{1} = [1,1]';  %<---use different initial point
y{1} = [2,3]';
y0{1} = [-1,1]';
g= @(x) A{1}(x)*y0{1}-y{1};
x=x0{1};

for m=1:10

Ja=@(x) [-4*sin(x(1)),     0;    %<----fix the calculus
         5*cos(x(1)),     0;
         -2*sin(x(1)+x(2)), -2*sin(x(1)+x(2));
          cos(x(1)+x(2)) ,  cos(x(1)+x(2))];

Jg=kron(y0{1}.',eye(2))*Ja(x);

x = x-Jg\g(x);  %<---dont use inv()

end

g(x)

ans = 2×1

1.0e-14 * 0.2665 -0.1332

[Rigon]

Thank you so much Matt, but the thing is the initial point was given in my exercise, I don't think I can use different one.

[Torsten]

Then you will have to live with the singular jacobian at the initial point.

[Matt J]

It will never work. You can see analytically that the Jacobain at x0=[0;0] is singular. Something to ask your instructor, I suppose...

Or you can try with pinv()

x = x-pinv(Jg)*g(x);  %<---dont use inv()

[Rigon]

Thank you so much once again. One more thing I tried my code while changing the initial point to [-1,1] for the first case A{1} and it actually works I got the same result for x just like fsolve. How can I make it general to work for other two cases, also not sure about the g(x) because where is called my the code is

x =nlinSolve(RootProblem,x0{i});

RootProblem = @(z) A{i}(z)*y0{i}-y{i}; , so the variable is z not x. Sorry for this stupid question?

function x = nlinSolve (g,x0)


 
select = @(g,r,c) g(r,c);
G1 = @(x)select(g(x), 1, 1);    
G2 =  @(x)select(g(x), 2, 1);   
G3 =  @(x)select(g(x), 3, 1);
G4 =  @(x)select(g(x), 4, 1);
G5 =  @(x)select(g(x), 5, 1);
x = x0;         
for m=1:10 
%calculating Jacobian
for k=1:numel(x)
    for n =1: numel(x)
        Jakobian(k,n)= Derivative(eval(['G',num2str(k)]),x,x(n),n);
    end
end
x = x-Jakobian\g(x);  
end
end

The code where it is supposed to be called my function is:

score = 0;
error = '';
details = struct();
try
    % Iteration ueber testcases
     for i = 1:3
        % Aufstellen des nullstellenproblems
        RootProblem = @(z) A{i}(z)*y0{i}-y{i};
        T = 0;
        % Mitteln ueber Zeit
		timer = tic();
        for j = 1:10
            x = nlinSolve(RootProblem,x0{i});  % <--- my function is called here
        end
        T = toc(timer)/10;
        output = A{i}(x)*y0{i}-y{i};
        Err = norm(output);
        Sc = (Err+0.01)*(T+0.5);
        score = score+Sc;
        
        details.(['time',num2str(i)]) = T;
        details.(['error',num2str(i)]) = Err;
     end
     score = 0.017 / score;
catch err
    score = 0;
    error = err.message;
end

[Matt J]

eval(['G',num2str(k)])

Bad practice! This should be G{1}, G{2}, ...G{n}.

details.(['time',num2str(i)]) = T;

details.(['error',num2str(i)]) = Err;

Also bad. This should be,

details.time(i)=T;

details.error(i)=Err;

or

details(i).time=T;

details(i).error=Err;

how can I make it general to work for other two cases

If you have the symbolic toolbox, you could set up the function and its Jacobian generically like in the following,

n=2; %dimension

syms x [n,1]
Asym =   [4*cos(x1),2*cos(x1+x2);  %symbolic function
         5*sin(x1),1*sin(x1+x2)]

Asym = 

Jsym=jacobian(Asym(:))  %symbolic Jacobian

Jsym = 

A=matlabFunction(Asym); %generate numerical implementations
Ja=matlabFunction(Jsym);

A=@(x) wrapIt(x,A);
Ja=@(x) wrapIt(x,Ja);

%test
x0=[pi;0];
A(x0)

ans = 2×2

-4.0000 -2.0000 0.0000 0.0000

Ja(x0)

ans = 4×2

-0.0000 0 -5.0000 0 -0.0000 -0.0000 -1.0000 -1.0000

function out=wrapIt(x,fun)

       x=num2cell(x);
       out=fun(x{:});
end

[Rigon]

Sorry Matt for probably being unclear. So this score function was written from the instructor and it is suppose to check automatically if the algorithm written from students is correct ( it is correct if the score is close to 1). So I only have to write a function which is going to solve a nonlinear system of equations. Using the same input in order for the Score function to work

function x = nlinSolve (g,x0)

Also I have 3 different cases, which you can see at the very beginning. So I am not suppose to change the Score4 function, just to write this nlinSove

[Rigon]

So I was thinking is it posible to define a function without using a variable here for example

G1 = @(x)select(g(x), 1, 1);

g(x) because g can be another function which is RootProblem = @(z) A{i}(z)*y0{i}-y{i}; using the variable z. So clearly I am doing that wrong. So how to select a value of function without using the variable x but something that works for whatever variable. I hope I didn't complicate this. Thank you because I am taking so much of your time.

[Matt J]

The way I see it, you have two choices:

(1) You approximate the Jacobian with finite differences.

(2) You apply symbolic math as discussed above. However, you would have to write a function, call it say func2sym(),that converts the anonymous functions passed to you by Score4 into a symbolic expression. There are ways you can do that using func2str() in conjunction with str2sym(), but it takes a little bit of thinking how to do that.

A = @(x)[4*cos(x(1)),2*cos(x(1)+x(2));
         5*sin(x(1)),1*sin(x(1)+x(2))]

A = function_handle with value:
    @(x)[4*cos(x(1)),2*cos(x(1)+x(2));5*sin(x(1)),1*sin(x(1)+x(2))]

Asym=func2sym(A)

Asym = 

whos A Asym

  Name      Size            Bytes  Class              Attributes

  A         1x1                32  function_handle              
  Asym      2x2                 8  sym                          

[Rigon]

Is there anyway I can make this code simplier but that still can work probably deleting some lines that are not that important? If not, I am okay with this one, still happy that it works.

Thank you once more Matt, you were very helpful really and I appreciate that.

[Matt J]

Thank you once more Matt, you were very helpful really and I appreciate that.

You are quite welcome, but if your question has been resolved, please Accept-click the answer.

Is there anyway I can make this code simplier but that still can work probably deleting some lines that are not that important?

I'm not sure you need the try...catch in funwrapper. For the purpose of this exercise, you know the function handles you are given will not be capable of returning their own Jacobian calculations. I'm also not sure why the Jacobian scaling has been introduced. Your first implementation was plain-vanilla Newton's method and it seemed to work fine.

[Rigon]

Matt you are a legend, keep up the good work.

[munazama]

1-cos3x1 + cos3x2 = 0

1-cos5x1 + cos5x2 = 0

can anybody plz provide its coding and solution using newton method. thanks in advance