Solution of nonlinear equation is different for different initial values

I have two codes. The first one is for two equations and the 2nd one is for 5 equations. the first one is giving the perfect result. Even if we change the initial values, it gives the actual result.

The first code

% Define the functions
f1 = @(x1, x2) x1 + x1.*x2 - 4;
f2 = @(x1, x2) x1 + x2 - 3;

% Define the partial derivatives
df1_dx1 = @(x1, x2) 1 + x2;
df1_dx2 = @(x1, x2) x1;
df2_dx1 = @(x1, x2) 1;
df2_dx2 = @(x1, x2) 1;

% Initialize the iteration counter and the initial values for x1 and x2
i = 0;
x0 = [0; 0.59]; % Column vector
tolerance = 1e-6;

% Display header for results
disp('Iteration | x1        | x2');

Iteration | x1        | x2

% Iterative process
while true
    % Compute the Jacobian matrix
    J = [df1_dx1(x0(1), x0(2)), df1_dx2(x0(1), x0(2));
         df2_dx1(x0(1), x0(2)), df2_dx2(x0(1), x0(2))];
    
    % Compute the function values vector
    F = [f1(x0(1), x0(2)); 
         f2(x0(1), x0(2))];
    
    % Compute the new values of x1 and x2
    x_new = x0 - J\F; % Equivalent to inv(J)*F but more efficient
    
    % Display the current iteration and the values
    fprintf('%9d | %f | %f\n', i, x_new(1), x_new(2));
    
    % Check if the differences are small enough to stop
    if abs(x_new - x0) < tolerance
        break;
    end
    
    % Update for the next iteration
    x0 = x_new;
    i = i + 1;
end

        0 | 2.515723 | 0.484277
        1 | 2.257862 | 0.742138
        2 | 2.128931 | 0.871069
        3 | 2.064465 | 0.935535
        4 | 2.032233 | 0.967767
        5 | 2.016116 | 0.983884
        6 | 2.008058 | 0.991942
        7 | 2.004029 | 0.995971
        8 | 2.002015 | 0.997985
        9 | 2.001007 | 0.998993
       10 | 2.000504 | 0.999496
       11 | 2.000252 | 0.999748
       12 | 2.000126 | 0.999874
       13 | 2.000063 | 0.999937
       14 | 2.000031 | 0.999969
       15 | 2.000016 | 0.999984
       16 | 2.000008 | 0.999992
       17 | 2.000004 | 0.999996
       18 | 2.000002 | 0.999998
       19 | 2.000001 | 0.999999

But result regarding the second one changes with initial values though I followed the same way to write both the codes. What are the problems with the second code? How should I correct it? I request your suggestions.

(For the second one I should get the result as [0.995675432 0.0699 0.279 0.279 1.057e-12] ).

Thanks in advance.

My 2nd code:

M_0 = 5;
gamma = 5.00 / 3.00;
p1 = -1.86 * 10^-1;
p2 = 1.86 * 10^-1;
p3 = 2.8 * 10^-1;
p4 = 0.0151 ; % in radians
p5 = 0.35;
theta_ap = acos(p4);
theta_am = pi - theta_ap;

m1 = 3/2 + 1/(M_0^2*(gamma -1));


% Define the functions
f1 = @(x1, x2, x3, x4, x5) (1 + cos(theta_ap) / M_0) * x4 + (1 - cos(theta_am) / M_0) * x5 + x2 - sin(x1) * x3 - (p5 + p1);
f2 = @(x1, x2, x3, x4, x5) (sin(theta_ap) / M_0) * x4 - (sin(theta_am) / M_0) * x5 + cos(x1) * x3 - p2;
f3 = @(x1, x2, x3, x4, x5) (1 + 2 * cos(theta_ap) / M_0 + 1 / M_0^2) * x4 + (1 - 2 * cos(theta_am) / M_0 + 1 / M_0^2) * x5 + x2 - ...
           2 * sin(x1) * x3 - (p5 + 2 * p1 + p3 / M_0^2);
f4 = @(x1, x2, x3, x4, x5) (m1 * cos(theta_ap) / M_0 + 1 / 2 + gamma / ((gamma - 1) * M_0^2)) * x4 + ...
           (1 / 2 + gamma / ((gamma - 1) * M_0^2) - m1 * cos(theta_am) / M_0) * x5 + ...
           1 / 2 * x2 - m1 * sin(x1) * x3 - (m1 * p1 + gamma / (gamma - 1) * p3 / M_0^2 + 1 / 2 * p5);
f5 = @(x1, x2, x3, x4, x5) ((x4 + x5) * cos(theta_ap) / M_0 - x3 * sin(x1) - p1);

% Define the partial derivatives
df1_dx1 = @(x1, x2, x3, x4, x5) -cos(x1) * x3;
df1_dx2 = @(x1, x2, x3, x4, x5) 1;
df1_dx3 = @(x1, x2, x3, x4, x5) -sin(x1);
df1_dx4 = @(x1, x2, x3, x4, x5) 1 + cos(theta_ap)/M_0;
df1_dx5 = @(x1, x2, x3, x4, x5) 1 - cos(theta_am)/M_0;

df2_dx1 = @(x1, x2, x3, x4, x5) -sin(x1) * x3;
df2_dx2 = @(x1, x2, x3, x4, x5) 0;
df2_dx3 = @(x1, x2, x3, x4, x5) cos(x1);
df2_dx4 = @(x1, x2, x3, x4, x5) sin(theta_ap)/M_0;
df2_dx5 = @(x1, x2, x3, x4, x5) -sin(theta_am)/M_0;

df3_dx1 = @(x1, x2, x3, x4, x5) -2*cos(x1)*x3;
df3_dx2 = @(x1, x2, x3, x4, x5) 1;
df3_dx3 = @(x1, x2, x3, x4, x5) -2*sin(x1);
df3_dx4 = @(x1, x2, x3, x4, x5) 1 + 2*cos(theta_ap)/M_0 + 1/M_0^2;
df3_dx5 = @(x1, x2, x3, x4, x5) 1 - 2*cos(theta_am)/M_0 + 1/M_0^2;

df4_dx1 = @(x1, x2, x3, x4, x5) -(3/2 + 1/(gamma - 1) * 1/M_0^2) * cos(x1)*x3;
df4_dx2 = @(x1, x2, x3, x4, x5) 1/2;
df4_dx3 = @(x1, x2, x3, x4, x5) -(3/2 + 1/(gamma - 1) * 1/M_0^2) * cos(x1);
df4_dx4 = @(x1, x2, x3, x4, x5) m1 * cos(theta_ap)/M_0 + 1/2 + gamma/((gamma - 1)*M_0^2);
df4_dx5 = @(x1, x2, x3, x4, x5) 1/2 + gamma/((gamma - 1)*M_0^2) - m1 * cos(theta_am)/M_0;

df5_dx1 = @(x1, x2, x3, x4, x5) -x3 * cos(x1);
df5_dx2 = @(x1, x2, x3, x4, x5) 0;
df5_dx3 = @(x1, x2, x3, x4, x5) -sin(x1);
df5_dx4 = @(x1, x2, x3, x4, x5) cos(theta_ap)/M_0;
df5_dx5 = @(x1, x2, x3, x4, x5) cos(theta_ap)/M_0;



eps = 1e-6; %tolerance

% Initialize the iteration counter and the initial values
i = 0;
x0 = [0.6981, 0.3, 0.4, 0.7, 0.1]; % initial values


% Display header for results
disp('Iteration | x1        | x2           | x3         | x4         | x5');

Iteration | x1        | x2           | x3         | x4         | x5

% Iterative process
while true
    % Compute the Jacobian matrix
     % Assuming dfn_dxm denotes the partial derivative of function n with respect to variable m
     J = [df1_dx1(x0(1), x0(2), x0(3), x0(4), x0(5)), df1_dx2(x0(1), x0(2), x0(3), x0(4), x0(5)), df1_dx3(x0(1), x0(2), x0(3), x0(4), x0(5)), df1_dx4(x0(1), x0(2), x0(3), x0(4), x0(5)), df1_dx5(x0(1), x0(2), x0(3), x0(4), x0(5));
          df2_dx1(x0(1), x0(2), x0(3), x0(4), x0(5)), df2_dx2(x0(1), x0(2), x0(3), x0(4), x0(5)), df2_dx3(x0(1), x0(2), x0(3), x0(4), x0(5)), df2_dx4(x0(1), x0(2), x0(3), x0(4), x0(5)), df2_dx5(x0(1), x0(2), x0(3), x0(4), x0(5));
          df3_dx1(x0(1), x0(2), x0(3), x0(4), x0(5)), df3_dx2(x0(1), x0(2), x0(3), x0(4), x0(5)), df3_dx3(x0(1), x0(2), x0(3), x0(4), x0(5)), df3_dx4(x0(1), x0(2), x0(3), x0(4), x0(5)), df3_dx5(x0(1), x0(2), x0(3), x0(4), x0(5));
          df4_dx1(x0(1), x0(2), x0(3), x0(4), x0(5)), df4_dx2(x0(1), x0(2), x0(3), x0(4), x0(5)), df4_dx3(x0(1), x0(2), x0(3), x0(4), x0(5)), df4_dx4(x0(1), x0(2), x0(3), x0(4), x0(5)), df4_dx5(x0(1), x0(2), x0(3), x0(4), x0(5));
          df5_dx1(x0(1), x0(2), x0(3), x0(4), x0(5)), df5_dx2(x0(1), x0(2), x0(3), x0(4), x0(5)), df5_dx3(x0(1), x0(2), x0(3), x0(4), x0(5)), df5_dx4(x0(1), x0(2), x0(3), x0(4), x0(5)), df5_dx5(x0(1), x0(2), x0(3), x0(4), x0(5))];


   
    % Compute the function values vector
   F =   [f1(x0(1), x0(2), x0(3), x0(4), x0(5)); 
          f2(x0(1), x0(2), x0(3), x0(4), x0(5)); 
          f3(x0(1), x0(2), x0(3), x0(4), x0(5)); 
          f4(x0(1), x0(2), x0(3), x0(4), x0(5)); 
          f5(x0(1), x0(2), x0(3), x0(4), x0(5))];

    % Compute the new values of x1 and x2
    % x_n = x0 - J\F; % Equivalent to inv(J)*F but more efficient
    x_n = x0 - inv(J)*F;
    
    % Display the current iteration and the values
    fprintf('%9d | %f | %f  | %f |  | %f |  | %f\n', i, x_n(1), x_n(2), x_n(3), x_n(4), x_n(5));
    
    % Check if the differences are small enough to stop
    if abs(x_n - x0) < eps
       break;
    end

    
    % Update for the next iteration
    x0 = x_n;
    i = i + 1
    % result should come like [0.995675432  0.0699 0.279 0.279 1.057e-12];
end

        0 | 0.468812 | 0.468100  | 0.698100 |  | -0.310394 |  | 1.186594

i = 1

        1 | 0.262359 | 0.070000  | 0.698100 |  | -1.114889 |  | 1.394889

i = 2

        2 | 0.270942 | 0.070000  | 0.698100 |  | -1.076783 |  | 1.356783

i = 3

        3 | 0.270952 | 0.070000  | 0.698100 |  | -1.076716 |  | 1.356716

i = 4

        4 | 0.270952 | 0.070000  | 0.698100 |  | -1.076716 |  | 1.356716

Answer 1

[Torsten]

Most probably, your second system of 5 equations has multiple solutions.

If you add the lines

x0 = x_n;
F =   [f1(x0(1), x0(2), x0(3), x0(4), x0(5)); 
    f2(x0(1), x0(2), x0(3), x0(4), x0(5)); 
    f3(x0(1), x0(2), x0(3), x0(4), x0(5)); 
    f4(x0(1), x0(2), x0(3), x0(4), x0(5)); 
    f5(x0(1), x0(2), x0(3), x0(4), x0(5))];
norm(F)

at the end of your second code, you will see that the obtained x-vector solves your system (although it's not the solution you expected to get).

This gives a solution with all its components >=0:

M_0 = 5;
gamma = 5.00 / 3.00;
p1 = -1.86 * 10^-1;
p2 = 1.86 * 10^-1;
p3 = 2.8 * 10^-1;
p4 = 0.0151 ; % in radians
p5 = 0.35;
theta_ap = acos(p4);
theta_am = pi - theta_ap;
m1 = 3/2 + 1/(M_0^2*(gamma -1));
% Define the functions
f1 = @(x1, x2, x3, x4, x5) (1 + cos(theta_ap) / M_0) * x4 + (1 - cos(theta_am) / M_0) * x5 + x2 - sin(x1) * x3 - (p5 + p1);
f2 = @(x1, x2, x3, x4, x5) (sin(theta_ap) / M_0) * x4 - (sin(theta_am) / M_0) * x5 + cos(x1) * x3 - p2;
f3 = @(x1, x2, x3, x4, x5) (1 + 2 * cos(theta_ap) / M_0 + 1 / M_0^2) * x4 + (1 - 2 * cos(theta_am) / M_0 + 1 / M_0^2) * x5 + x2 - ...
           2 * sin(x1) * x3 - (p5 + 2 * p1 + p3 / M_0^2);
f4 = @(x1, x2, x3, x4, x5) (m1 * cos(theta_ap) / M_0 + 1 / 2 + gamma / ((gamma - 1) * M_0^2)) * x4 + ...
           (1 / 2 + gamma / ((gamma - 1) * M_0^2) - m1 * cos(theta_am) / M_0) * x5 + ...
           1 / 2 * x2 - m1 * sin(x1) * x3 - (m1 * p1 + gamma / (gamma - 1) * p3 / M_0^2 + 1 / 2 * p5);
f5 = @(x1, x2, x3, x4, x5) ((x4 + x5) * cos(theta_ap) / M_0 - x3 * sin(x1) - p1);
f = @(x)[f1(x(1),x(2),x(3),x(4),x(5));f2(x(1),x(2),x(3),x(4),x(5));f3(x(1),x(2),x(3),x(4),x(5));f4(x(1),x(2),x(3),x(4),x(5));f5(x(1),x(2),x(3),x(4),x(5))];
x0 = [0.995675432  0.0699 0.279 0.279 1.057e-12];
x = lsqnonlin(f,x0,zeros(5,1),inf(5,1),optimset('TolFun',1e-12,'TolX',1e-12))

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

x = 1×5

0.9240 0.0700 0.2341 0.2523 0.0277

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

norm(f(x))

ans = 4.1161e-14