Solution of nonlinear equation is different for different initial values
1 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(1 个)
Torsten
2024-3-29
编辑:Torsten
2024-3-29
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 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Gamma Functions 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!