Solving equation with two for loops
4 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
DS
am 15 Nov. 2022
Beantwortet: Walter Roberson
am 18 Nov. 2022
Hi I would like to find the values of gf for different r and H values. It is obtained by solving the equation S1 which includes gf on both left and right sides. Could you please help me in solving this?
clear all;
close all;
iiy = 1.8;
rho_s = 2160;
rho_w = 997;
M_s = 58.44 * 10^-3;
M_w = 18* 10^-3
o = 0.072;
H = linspace(0.1,1,100)
T = 22+273;
Rv = 8.314;
r = 1E-9:100E-9:200E-9;
syms gf
for i = 1:length(H)
for j =1:length(r)
S1 = ( gf -((rho_s/rho_w) *(1 + (M_w/M_s * iiy * H(i)./(( exp(2*o)./(rho_w*Rv*T.*gf*r(j))) - H(i)))))^(1/3));
S(j,i) = vpasolve(S1(j,i)==0,gf, [0 Inf] );
end
end
0 Kommentare
Akzeptierte Antwort
Walter Roberson
am 18 Nov. 2022
I suggest you switch to numeric solutions, and to providing solution hints. Some of the tests I did hinted there might be multiple solutions.
The bumps along the yellow ridge appear to be inherent somehow -- if you zoom in you will see "ribs" on the plot.
iiy = 1.8;
rho_s = 2160;
rho_w = 997;
M_s = 58.44 * 10^-3;
M_w = 18* 10^-3 ;
M_w = 0.0180;
o = 0.072;
H = linspace(0.1,1,100);
T = 22+273;
Rv = 8.314;
r = linspace(1E-9,1000E-9, 50);
numH = length(H);
numr = length(r);
S = NaN(numr, numH);
guess = 1.5;
opt = optimoptions('fsolve','Display','off');
for i = 1:numH
for j = 1:numr
S1 = @(gf) ( gf -((rho_s/rho_w) *(1 + (M_w/M_s * iiy * H(i)./(( exp(2*o)./(rho_w*Rv*T.*gf*r(j))) - H(i)))))^(1/3));
SS = fsolve(S1, guess, opt );
if ~isempty(SS) && imag(SS) == 0
S(j,i) = double(SS);
guess = SS;
end
end
end
surf(H, r, S, 'edgecolor', 'none')
xlabel('H'); ylabel('r'); zlabel('S')
view([57,39])
0 Kommentare
Weitere Antworten (2)
Voss
am 15 Nov. 2022
S1 is calculated for a given i and j, so it is not correct to index it with i and j on the next line.
That is, use S1==0 instead of S1(j,i) == 0
clear all;
close all;
iiy = 1.8;
rho_s = 2160;
rho_w = 997;
M_s = 58.44 * 10^-3;
M_w = 18* 10^-3
o = 0.072;
H = linspace(0.1,1,100)
T = 22+273;
Rv = 8.314;
r = 1E-9:100E-9:200E-9;
syms gf
for i = 1:length(H)
for j =1:length(r)
S1 = ( gf -((rho_s/rho_w) *(1 + (M_w/M_s * iiy * H(i)./(( exp(2*o)./(rho_w*Rv*T.*gf*r(j))) - H(i)))))^(1/3));
% S(j,i) = vpasolve(S1(j,i)==0,gf, [0 Inf] );
S(j,i) = vpasolve(S1==0,gf, [0 Inf] );
end
end
disp(S)
3 Kommentare
Walter Roberson
am 15 Nov. 2022
The implication is that there is a combination of values that vpasolve is not able to find a solution for.
Walter Roberson
am 15 Nov. 2022
iiy = 1.8;
rho_s = 2160;
rho_w = 997;
M_s = 58.44 * 10^-3;
M_w = 18* 10^-3
o = 0.072;
H = linspace(0.1,1,100)
T = 22+273;
Rv = 8.314;
r = 1E-9:100E-9:200E-9;
syms gf
numH = length(H);
numr = length(r);
S = zeros(numr, numH);
for i = 1:numH
for j = 1:numr
S1 = ( gf -((rho_s/rho_w) *(1 + (M_w/M_s * iiy * H(i)./(( exp(2*o)./(rho_w*Rv*T.*gf*r(j))) - H(i)))))^(1/3));
S(j,i) = vpasolve(S1==0, gf, [0 Inf] );
end
end
format long g
S
3 Kommentare
Torsten
am 15 Nov. 2022
Bearbeitet: Torsten
am 15 Nov. 2022
Your equation might have no solution or vpasolve might fail to compute it ...
iiy = 1.8;
rho_s = 2160;
rho_w = 997;
M_s = 58.44 * 10^-3;
M_w = 18* 10^-3
M_w = 0.0180
o = 0.072;
H = linspace(0.1,1,100)
T = 22+273;
Rv = 8.314;
r = linspace(1E-9,1000E-9, 20);
syms gf
numH = length(H);
numr = length(r);
S = NaN(numr, numH);
for i = 1:numH
for j = 1:numr
S1 = ( gf -((rho_s/rho_w) *(1 + (M_w/M_s * iiy * H(i)./(( exp(2*o)./(rho_w*Rv*T.*gf*r(j))) - H(i)))))^(1/3));
SS = vpasolve(S1==0, gf, [0 Inf] );
if ~isempty(SS)
S(j,i) = double(SS);
end
end
end
format long g
S
Siehe auch
Kategorien
Mehr zu Loops and Conditional Statements 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!