Solving equation with two for loops

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 

M_w = 0.0180

o = 0.072;
H = linspace(0.1,1,100)

H = 1×100

0.1000 0.1091 0.1182 0.1273 0.1364 0.1455 0.1545 0.1636 0.1727 0.1818 0.1909 0.2000 0.2091 0.2182 0.2273 0.2364 0.2455 0.2545 0.2636 0.2727 0.2818 0.2909 0.3000 0.3091 0.3182 0.3273 0.3364 0.3455 0.3545 0.3636

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

Error using sub2ind
Out of range subscript.

Error in indexing (line 1074)
R_tilde = sub2ind(size(L), Idx.subs{:});

Answer 1

[Walter Roberson]

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])

Answer 2

[Voss]

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 

M_w = 0.0180

o = 0.072;
H = linspace(0.1,1,100)

H = 1×100

0.1000 0.1091 0.1182 0.1273 0.1364 0.1455 0.1545 0.1636 0.1727 0.1818 0.1909 0.2000 0.2091 0.2182 0.2273 0.2364 0.2455 0.2545 0.2636 0.2727 0.2818 0.2909 0.3000 0.3091 0.3182 0.3273 0.3364 0.3455 0.3545 0.3636

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)

[DS]

r = linspace(1E-9,1000E-9, 200);

If I use this one for the value of r,

why am I getting this error?

Unable to perform assignment because the size of the left side is 1-by-1 and the size of

the right side is 0-by-1.

Error in sym/privsubsasgn (line 1229)

L_tilde2 = builtin('subsasgn',L_tilde,struct('type','()','subs',{varargin}),R_tilde);

Error in sym/subsasgn (line 1060)

C = privsubsasgn(L,R,inds{:});

Error in lietal (line 34)

S(j,i) = vpasolve(S1==0,gf, [0 Inf] );

>>

[Walter Roberson]

The implication is that there is a combination of values that vpasolve is not able to find a solution for.

[DS]

How can I solve this? Should I use some other functions?

Answer 3

[Walter Roberson]

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)

H = 1×100

0.1000 0.1091 0.1182 0.1273 0.1364 0.1455 0.1545 0.1636 0.1727 0.1818 0.1909 0.2000 0.2091 0.2182 0.2273 0.2364 0.2455 0.2545 0.2636 0.2727 0.2818 0.2909 0.3000 0.3091 0.3182 0.3273 0.3364 0.3455 0.3545 0.3636

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

S = 2×100

1.29402159925611 1.29402755858776 1.2940335182164 1.29403947814207 1.29404543836478 1.29405139888457 1.29405735970146 1.29406332081548 1.29406928222665 1.294075243935 1.29408120594055 1.29408716824334 1.29409313084338 1.29409909374071 1.29410505693534 1.29411102042732 1.29411698421665 1.29412294830337 1.2941289126875 1.29413487736908 1.29414084234812 1.29414680762465 1.2941527731987 1.2941587390703 1.29416470523946 1.29417067170622 1.29417663847061 1.29418260553264 1.29418857289235 1.29419454054976 1.30076236019233 1.30140057925328 1.3020421735055 1.30268717464784 1.30333561482099 1.30398752661578 1.30464294308161 1.30530189773512 1.3059644245691 1.30663055806151 1.30730033318485 1.30797378541565 1.3086509507442 21.4298073297381 1.31001656728512 1.31070509313834 1.31139748139228 1.31209377076126 1.31279400053724 1.31349821060138 1.31420644143594 1.31491873413649 1.31563513042437 1.31635567265952 1.31708040385366 1.31780936768376 1.31854260850588 1.31928017136941 1.32002210203163 1.32076844697273

[DS]

iiy = 1.8;
rho_s = 2160;
rho_w = 997; 
M_s = 58.44 * 10^-3; 
M_w = 18* 10^-3 

M_w = 0.0180

M_w = 0.0180

M_w = 0.0180

o = 0.072;
H = linspace(0.1,1,100)

H = 1×100

0.1000 0.1091 0.1182 0.1273 0.1364 0.1455 0.1545 0.1636 0.1727 0.1818 0.1909 0.2000 0.2091 0.2182 0.2273 0.2364 0.2455 0.2545 0.2636 0.2727 0.2818 0.2909 0.3000 0.3091 0.3182 0.3273 0.3364 0.3455 0.3545 0.3636

T = 22+273;
Rv = 8.314;
r = linspace(1E-9,1000E-9, 20);
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

Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 0-by-1.

format long g
S

[DS]

Unable to perform assignment because the size of the left side is 1-by-1 and the size of

the right side is 0-by-1.

I still got this error. Could you please help me what should I do?

[Torsten]

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