Solution of nonlinear equation is different for different initial values
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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');
% 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
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');
% 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 Kommentare
Antworten (1)
Torsten
am 29 Mär. 2024
Bearbeitet: Torsten
am 29 Mär. 2024
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))
norm(f(x))
0 Kommentare
Siehe auch
Kategorien
Mehr zu Gamma Functions finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!