finding non zero solution for homogeneous system of equations

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');

Answer 1

[Aquatris]

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),:)

nonzeroSolutions = 557x7

-11.1439 -2.4841 0.5966 -0.3815 0.1024 -0.1500 0.0153 -1.1769 1.1402 0.9476 -0.1999 -0.1267 0.0423 -0.1175 -161.9110 -3.7509 -0.9152 15.7362 3.5740 -39.5678 39.4107 -5.9693 -4.2321 0.1918 0.0064 -0.0264 -0.8981 -0.0094 -6.0161 0.2271 -0.9818 -4.0342 -0.5380 1.5874 -1.7717 -7.1687 5.3215 2.1127 -0.2700 -0.7532 0.3789 -0.1887 -10.5027 -7.6824 -1.8260 -0.1935 -1.0856 -0.5672 0.1584 -11.6710 -8.5242 1.4140 0.1407 -1.1966 -0.6349 -0.1218 -11.7617 0.0084 10.3163 -0.0083 0.4097 -0.0118 2.9077 -12.3794 0.2234 10.8350 0.0211 0.2714 0.0482 2.8422

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

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

%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');

[John D'Errico]

+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.