finding non zero solution for homogeneous system of equations
4 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
The following code is giving constant values of Nusselt number. The solutions are approximately close to zero.
D = 0.1;
L = 0.1;
B = 0.1;
xi = 0.1;
sigma1 = B^3 .* D^2 .* L^3;
sigma2 = D^2 .* L^3;
Ra_values = (10:100:100000).';
solutions = zeros(length(Ra_values), 7);
%Nu = zeros(length(Ra_values),1);
for i = 1:length(Ra_values)
Ra = Ra_values(i);
R = Ra * xi;
f1 = @(A) -((5.*A(2).*A(6).*pi^5.*B^2.*D^2.*L^4)./8 + (A(2).*A(6).*pi^5.*D^4.*L^4)./4)./sigma1 - (B .* ((A(3).*A(7).*D^4.*pi^5)/4 - (A(1).*L^4.*pi^4)./2 - (32.*A(5).*L^4.*Ra)./9 + (A(1).*L^4.*R.*pi^2)./2 + (5.*A(3).*A(7).*D^2.*L^2.*pi^5)./8))./sigma2;
f2 = @(A) -(B^2 .* (A(2).*A(5).*D^2.*L^4.*pi^5 - (A(2).*D^2.*L^4.*pi^4)./2 - (A(1).*A(6).*D^2.*L^4.*pi^5)./8 - (4.*A(6).*D^2.*L^4.*Ra)./9 + (A(4).*A(7).*D^2.*L^4.*pi^5)./2 + (3.*A(4).*A(7).*D^4.*L^2.*pi^5)./16 + (A(2).*D^2.*L^4.*R.*pi^2)./4) - (A(2).*D^4.*L^4.*pi^4)./4)./sigma1 - (B .* ((A(4).*A(7).*D^4.*pi^5)./8 - (A(2).*L^4.*pi^4)./4 - (16.*A(6).*L^4.*Ra)./9 + (A(2).*L^4.*R.*pi^2)./4 + (5.*A(4).*A(7).*D^2.*L^2.*pi^5)./16))./sigma2;
f3 = @(A) (B .* ((16.*A(7).*L^4.*Ra)./9 + (A(3).*D^4.*pi^4)./4 + (A(3).*L^4.*pi^4)./4 - (A(3).*L^4.*R.*pi^2)./4 + (A(3).*D^2.*L^2.*pi^4)./2 + (A(1).*A(7).*D^2.*L^2.*pi^5)./8 - A(3).*A(5).*D^2.*L^2.*pi^5 - (A(4).*A(6).*D^2.*L^2.*pi^5)./2))./sigma2 - (B^2 .* ((3.*A(4).*A(6).*pi^5.*D^4.*L^2)./16 + (5.*A(4).*A(6).*pi^5.*D^2.*L^4)./16) + (A(4).*A(6).*D^4.*L^4.*pi^5)./8)./sigma1;
f4 = @(A) (B .* pi^2 .* (2.*A(4).*D^4.*pi^2 - 2.*A(4).*L^4.*R + 2.*A(4).*L^4.*pi^2 + 4.*A(4).*D^2.*L^2.*pi^2 + A(2).*A(7).*D^2.*L^2.*pi^3 - 8.*A(3).*A(6).*D^2.*L^2.*pi^3 - 8.*A(4).*A(5).*D^2.*L^2.*pi^3))./(16.*D^2.*L^3) - ((B^2 .* pi^2 .* (2.*A(4).*D^2.*L^4.*R - 4.*A(4).*D^2.*L^4.*pi^2 - 4.*A(4).*D^4.*L^2.*pi^2 + 8.*A(2).*A(7).*D^2.*L^4.*pi^3 + A(2).*A(7).*D^4.*L^2.*pi^3 - A(3).*A(6).*D^2.*L^4.*pi^3 + A(3).*A(6).*D^4.*L^2.*pi^3 + 8.*A(4).*A(5).*D^2.*L^4.*pi^3))./16 - (A(4).*D^4.*L^4.*pi^4)./8)./sigma1;
f5 = @(A) (B^2 .* ((pi^5.*A(2)^2.*D^2.*L^4)./4 + (pi^5.*A(4)^2.*D^4.*L^2)./4 + (pi^5.*A(4)^2.*D^2.*L^4)./8) + (A(2)^2.*D^4.*L^4.*pi^5)./4 + (A(4)^2.*D^4.*L^4.*pi^5)./8)./sigma1 + (B .* ((A(3)^2.*D^4.*pi^5)/4 + (A(4)^2.*D^4.*pi^5)./8 + (8.*A(1).*L^4.*Ra)./9 + 8.*A(5).*L^4.*pi^4 - 2.*A(5).*L^4.*R.*pi^2 + (A(3)^2.*D^2.*L^2.*pi^5)./4 + (A(4)^2.*D^2.*L^2.*pi^5)./8))./sigma2;
f6 = @(A) (B^2 .* ((4.*A(2).*D^2.*L^4.*Ra)./9 + 2.*A(6).*D^2.*L^4.*pi^4 + (A(1).*A(2).*D^2.*L^4.*pi^5)./8 + (A(3).*A(4).*D^2.*L^4.*pi^5)./16 + (3.*A(3).*A(4).*D^4.*L^2.*pi^5)./16 - (A(6).*D^2.*L^4.*R.*pi^2)./4) + (A(6).*D^4.*L^4.*pi^4)./4)./sigma1 + (B .* (4.*A(2).*L^4.*Ra)./9 + 4.*A(6).*L^4.*pi^4 + (A(3).*A(4).*D^4.*pi^5)./4 - A(6).*L^4.*R.*pi^2 + (A(3).*A(4).*D^2.*L^2.*pi^5)./4)./sigma2;
f7 = @(A) (B .* ((4.*A(3).*L^4.*Ra)./9 + (A(7).*D^4.*pi^4)./4 + 4.*A(7).*L^4.*pi^4 - A(7).*L^4.*R.*pi^2 + 2.*A(7).*D^2.*L^2.*pi^4 + (A(1).*A(3).*D^2.*L^2.*pi^5)./8 + (A(2).*A(4).*D^2.*L^2.*pi^5)./16))./sigma2 + (B^2 .* ((3.*A(2).*A(4).*pi^5.*D^4.*L^2)./16 + (A(2).*A(4).*pi^5.*D^2.*L^4)./4) + (A(2).*A(4).*D^4.*L^4.*pi^5)./4)./sigma1;
F = @(A) [f1(A); f2(A); f3(A); f4(A); f5(A); f6(A); f7(A)];
A0 = [0.01 0 0 0 0.01 0 0];
A = fsolve(F, A0);
solutions(i, :) = A;
F(A)
end
A1_1_1 = solutions(:, 1);
A2_1_1 = solutions(:, 5);
Nu = xi - 1./2 - (pi^3 ./ Ra_values) .* A1_1_1 - (8 .* pi^3 ./ Ra_values) .* A2_1_1;
%Nu(i) = xi - 1./2 - (pi^3 ./ Ra_values(i)) .* A(1) - (8 * pi^3 ./ Ra_values(i)) .* A(5);
%end
plot(Ra_values, Nu);
xlabel('Ra');
ylabel('Nu');
0 comentarios
Respuestas (1)
Aquatris
el 24 de Jun. de 2024
Editada: Aquatris
el 24 de Jun. de 2024
Since zeros(7,1) is kind of a solution for you, you should select your initial value away from that point. Below you can find a simple adjustment to your code that shows some of the nonzero solutions. If you want to find nonzero solutions for all Ra_values, then you should implement a simple algorithm that changes the initial values for specific Ra until a nonzero solution is found
D = 0.1;
L = 0.1;
B = 0.1;
xi = 0.1;
sigma1 = B^3 .* D^2 .* L^3;
sigma2 = D^2 .* L^3;
Ra_values = (10:100:100000).';
solutions = zeros(length(Ra_values), 7);
%Nu = zeros(length(Ra_values),1);
options = optimoptions('fsolve','Display','none');
for i = 1:length(Ra_values)
Ra = Ra_values(i);
R = Ra * xi;
f1 = @(A) -((5.*A(2).*A(6).*pi^5.*B^2.*D^2.*L^4)./8 + (A(2).*A(6).*pi^5.*D^4.*L^4)./4)./sigma1 - (B .* ((A(3).*A(7).*D^4.*pi^5)/4 - (A(1).*L^4.*pi^4)./2 - (32.*A(5).*L^4.*Ra)./9 + (A(1).*L^4.*R.*pi^2)./2 + (5.*A(3).*A(7).*D^2.*L^2.*pi^5)./8))./sigma2;
f2 = @(A) -(B^2 .* (A(2).*A(5).*D^2.*L^4.*pi^5 - (A(2).*D^2.*L^4.*pi^4)./2 - (A(1).*A(6).*D^2.*L^4.*pi^5)./8 - (4.*A(6).*D^2.*L^4.*Ra)./9 + (A(4).*A(7).*D^2.*L^4.*pi^5)./2 + (3.*A(4).*A(7).*D^4.*L^2.*pi^5)./16 + (A(2).*D^2.*L^4.*R.*pi^2)./4) - (A(2).*D^4.*L^4.*pi^4)./4)./sigma1 - (B .* ((A(4).*A(7).*D^4.*pi^5)./8 - (A(2).*L^4.*pi^4)./4 - (16.*A(6).*L^4.*Ra)./9 + (A(2).*L^4.*R.*pi^2)./4 + (5.*A(4).*A(7).*D^2.*L^2.*pi^5)./16))./sigma2;
f3 = @(A) (B .* ((16.*A(7).*L^4.*Ra)./9 + (A(3).*D^4.*pi^4)./4 + (A(3).*L^4.*pi^4)./4 - (A(3).*L^4.*R.*pi^2)./4 + (A(3).*D^2.*L^2.*pi^4)./2 + (A(1).*A(7).*D^2.*L^2.*pi^5)./8 - A(3).*A(5).*D^2.*L^2.*pi^5 - (A(4).*A(6).*D^2.*L^2.*pi^5)./2))./sigma2 - (B^2 .* ((3.*A(4).*A(6).*pi^5.*D^4.*L^2)./16 + (5.*A(4).*A(6).*pi^5.*D^2.*L^4)./16) + (A(4).*A(6).*D^4.*L^4.*pi^5)./8)./sigma1;
f4 = @(A) (B .* pi^2 .* (2.*A(4).*D^4.*pi^2 - 2.*A(4).*L^4.*R + 2.*A(4).*L^4.*pi^2 + 4.*A(4).*D^2.*L^2.*pi^2 + A(2).*A(7).*D^2.*L^2.*pi^3 - 8.*A(3).*A(6).*D^2.*L^2.*pi^3 - 8.*A(4).*A(5).*D^2.*L^2.*pi^3))./(16.*D^2.*L^3) - ((B^2 .* pi^2 .* (2.*A(4).*D^2.*L^4.*R - 4.*A(4).*D^2.*L^4.*pi^2 - 4.*A(4).*D^4.*L^2.*pi^2 + 8.*A(2).*A(7).*D^2.*L^4.*pi^3 + A(2).*A(7).*D^4.*L^2.*pi^3 - A(3).*A(6).*D^2.*L^4.*pi^3 + A(3).*A(6).*D^4.*L^2.*pi^3 + 8.*A(4).*A(5).*D^2.*L^4.*pi^3))./16 - (A(4).*D^4.*L^4.*pi^4)./8)./sigma1;
f5 = @(A) (B^2 .* ((pi^5.*A(2)^2.*D^2.*L^4)./4 + (pi^5.*A(4)^2.*D^4.*L^2)./4 + (pi^5.*A(4)^2.*D^2.*L^4)./8) + (A(2)^2.*D^4.*L^4.*pi^5)./4 + (A(4)^2.*D^4.*L^4.*pi^5)./8)./sigma1 + (B .* ((A(3)^2.*D^4.*pi^5)/4 + (A(4)^2.*D^4.*pi^5)./8 + (8.*A(1).*L^4.*Ra)./9 + 8.*A(5).*L^4.*pi^4 - 2.*A(5).*L^4.*R.*pi^2 + (A(3)^2.*D^2.*L^2.*pi^5)./4 + (A(4)^2.*D^2.*L^2.*pi^5)./8))./sigma2;
f6 = @(A) (B^2 .* ((4.*A(2).*D^2.*L^4.*Ra)./9 + 2.*A(6).*D^2.*L^4.*pi^4 + (A(1).*A(2).*D^2.*L^4.*pi^5)./8 + (A(3).*A(4).*D^2.*L^4.*pi^5)./16 + (3.*A(3).*A(4).*D^4.*L^2.*pi^5)./16 - (A(6).*D^2.*L^4.*R.*pi^2)./4) + (A(6).*D^4.*L^4.*pi^4)./4)./sigma1 + (B .* (4.*A(2).*L^4.*Ra)./9 + 4.*A(6).*L^4.*pi^4 + (A(3).*A(4).*D^4.*pi^5)./4 - A(6).*L^4.*R.*pi^2 + (A(3).*A(4).*D^2.*L^2.*pi^5)./4)./sigma2;
f7 = @(A) (B .* ((4.*A(3).*L^4.*Ra)./9 + (A(7).*D^4.*pi^4)./4 + 4.*A(7).*L^4.*pi^4 - A(7).*L^4.*R.*pi^2 + 2.*A(7).*D^2.*L^2.*pi^4 + (A(1).*A(3).*D^2.*L^2.*pi^5)./8 + (A(2).*A(4).*D^2.*L^2.*pi^5)./16))./sigma2 + (B^2 .* ((3.*A(2).*A(4).*pi^5.*D^4.*L^2)./16 + (A(2).*A(4).*pi^5.*D^2.*L^4)./4) + (A(2).*A(4).*D^4.*L^4.*pi^5)./4)./sigma1;
F = @(A) [f1(A); f2(A); f3(A); f4(A); f5(A); f6(A); f7(A)];
A0 = randn(7,1)*1e3;
A = fsolve(F, A0,options);
solutions(i, :) = A;
end
A1_1_1 = solutions(:, 1);
A2_1_1 = solutions(:, 5);
Nu = xi - 1./2 - (pi^3 ./ Ra_values) .* A1_1_1 - (8 .* pi^3 ./ Ra_values) .* A2_1_1;
nonzeroSolutions = solutions(all(abs(solutions)>1e-8,2),:)
%Nu(i) = xi - 1./2 - (pi^3 ./ Ra_values(i)) .* A(1) - (8 * pi^3 ./ Ra_values(i)) .* A(5);
%end
plot(Ra_values, Nu);
xlabel('Ra');
ylabel('Nu');
1 comentario
John D'Errico
el 24 de Jun. de 2024
Editada: John D'Errico
el 24 de Jun. de 2024
+1 of course. I very much appreciate seeing posts like yours. Well reasoned, good explanations and in good depth. I hope we continue to see you posting for a long time on the forum.
I would only suggest an idea that when you start with a random point at each pass for subtly different values of the parameter Ra, especially since this is a nonlinear system, you may be finding totally different solutions. It may be better to use the previous solution as a start point for fsolve at the next value of Ra.
Ver también
Categorías
Más información sobre Systems of Nonlinear Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!