MATLAB: Solve: Working only sometimes

MATLABsolve

Hello all,

I am using solve to find the numerical value of two constants from two polynomial equations. As I vary the number of terms in the equation (set by changing the value of nmax), solve will return a solution for nmax <= 35, but not for higher values, instead stating: Warning: Could not extract individual solutions. Returning a MuPAD set object.. Do you have any suggestions on how I might be able to find values for my constants at higher polynomial orders?

Many thanks!

More info: I have analytically solved a differential equation, and the solution is the sum of two power series with two unknowns. For certain values of my constants, I can truncate the series after only a few terms and get a convergent solution. However, as I vary the value of my constants, I need to go to higher orders to get a convergent solution.

Code:

% Specify the appropriate constants for finding the bounds of the depletion region

k = 8.6173324E-5; % The Bolzmann constant (eV/K)

T = 300; % The temperature (K)

q = 1.60219E-19; % The fundamental electrical charge (C)

VT = k*T; % The thermal voltage (eV)

NA = 5E16; % The p-type, core doping (cm^-3)

ND = 9E17; % The n-type, shell doping (cm^-3)

NC = 5.758E17; % Effective density of states in the conduction band (cm^-3)

NV = 9.543E18; % Effective density of states in the valence band (cm^-3)

EG = 1.6071; % The material bandgap (eV)

ni = sqrt(NC*NV)*exp(-EG/(2*VT)); % The intrinsic carrier density (cm^-3)

Phi_0 = VT*log(NA*ND/ni^2); % The built in voltage (eV)

V = 0; % The applied voltage (V)

d1 = 0.1E-4; % Emitter thickness (cm)

d2 = 1.9E-4; % The core radius (cm)

eps_GaP = 11.1; % The dielectric constant of GaP

eps_GaAs = 12.9; % The dielectric constant of GaAs

eps0 = 8.854187817620E-14; % Vacuum permittivty in F/cm

eps = (11.1 + 1.8*0.9)*eps0; % The overall permittivity

% Solve for the limits of the depletion region

syms x2 x4;

S = solve(Phi_0 + V == q*NA/(6*eps)*d2^2+q/(3*eps*d2)*(NA*x4^3+ND*x2^3)+q /(2*eps)*(ND*x2^2-NA*x4^2), NA*(d2^3-x4^3) == ND*((d2+x2)^3-d2^3));

for a = 1:length(S.x2)

    if isreal(S.x2(a)) && S.x2(a)>0 && isreal(S.x4(a)) && S.x4(a)>0
        x2_val = double(S.x2(a));
        x4_val = double(S.x4(a));
    end

end

% Specify the appropriate constants for finding the currents in the quasi-neutral regions

nmax = 40; % Set the sum limit to a finite, maximum value

Ln = 1E-4; % The diffusion length in the n-type region

Lp = 1E-6; % The diffusion length in the p-type region

mun = 7661; % The electron mobility (cm^2/V/sec)

mup = 367.5; % The hole mobility (cm^2/V/sec)

Dn = mun*k*T; % The diffusion coefficient for electrons (cm^2/sec)

Dp = mup*k*T; % The diffusion coefficient for holes (cm^2/sec)

N_im = 0.1; % The imaginary part of the index of refraction

wl = 6E-5; % The wavelength, in cm.

alpha = 4*pi*N_im/wl; % The absorption coefficent (cm^-1)

inc_pow = 100; % AM1.5G power (mW/cm^2)

h = 6.626E-34; % Planck's constant (Jsec)

c = 3E10 ; % Speed of light (cm/s)

flux = inc_pow/1000*wl/(c*h); % The incident photon flux (photons/cm^2/sec)

% Solve for the current in the quasi-neutral shell

bvals = sym('bvals',[nmax 1]); % Define a matrix to hold the symbolic values for the polynomial coefficients

syms b0 B; % Define the constants to be solved for with the boundary conditions

% Populate the coefficient matrix

bvals(1,1) = b0;

bvals(2,1) = -alpha*b0;

bvals(3,1) = -1/6*(4*alpha*bvals(2,1)+(alpha^2-1/Lp^2)*bvals(1,1)+alpha*flux/Dp);

for i = 4:nmax

    bvals(i,1) = 1/(i-1)/i*(2*alpha*(i-1)*bvals(i-1,1)+(alpha^2-1/Lp^2)*bvals(i-2,1));

end

% Calculate the various components of the minority carrier density

C2 = 0;

C3 = 0;

for k = 1:nmax

    C2 = (d2+x2_val)^(k-1)*bvals(k,1) + C2;
    C3 = (d1+d2)^(k-1)*bvals(k,1) + C3;

end

Hg2 = 0;

Hg3 = 0;

for l = 1:nmax

    Hg2 = (1/(Lp^2*(l+1)*(l+2)))^(l-1)*(d2+x2_val)^(2*(l-1)) + Hg2;
    Hg3 = (1/(Lp^2*(l+1)*(l+2)))^(l-1)*(d1+d2)^(2*(l-1)) + Hg3;

end

Sp = 100;

C2 = C2*exp(alpha*(x2_val-d1)) + B*Hg2;

C3 = C3 + B*Hg3;

dC3 = 0;

for p=1:nmax

    dCa = (p-1)*bvals(p,1)*(d2+d1)^(p-2);
    dCb = alpha*bvals(p,1)*(d2+d1)^(p-1);
    dCc = B*2*(p-1)*(1/(Lp^2*(p+1)*(p+2)))^(p-1)*(d1+d2)^(2*p-3);
    dC3 = dCa+dCb+dCc+dC3;

end

S3 = solve(C2 == ni^2/ND*(exp(V/VT)-1),-Dp*dC3==Sp*C3);

Best Answer

Some of your variable get 'Inf' in between. If you want to use variable precision arithmetic, make sure you use syms. I get a result if I change the following lines:

    Hg2 = sym(1/(Lp^2*(l+1)*(l+2)))^(l-1)*(d2+x2_val)^(2*(l-1)) + Hg2;
    Hg3 = sym(1/(Lp^2*(l+1)*(l+2)))^(l-1)*(d1+d2)^(2*(l-1)) + Hg3;

and

    dCc = B*2*(p-1)*sym(1/(Lp^2*(p+1)*(p+2)))^(p-1)*(d1+d2)^(2*p-3);

Not sure however if the result is correct:

 >> double(S3.B)
 ans =
   -7.0781e-07
 >> double(S3.b0)
 ans =
   6.4681e+08

Related Solutions

MATLAB: How to write a code for iteration

After the first box, you write

deltaT = epsilon + 1; % so that the while loop is entered
while deltaT >= epsilon
    % put code of the second box here
    deltaT = % add equation for new updated value of deltaT
 end

MATLAB: My fluid simulation exploding to very large numbers when I take small step values

I found a solution in the comments. I updated the value of n at each stage of the computations of j,k,l,o and now it is working as expected. Here is a code solution for the fluid simulation that works as intended:

clear;
clc;
format long g;
%% Defining the constants
L = 0.6; % Characteristic Length of evaporation surface in m
W = 0.6; % Width of the pan in m
A = L*W; % Surface area of one pan in m^2
d = 0.032; % Height of air columns in m 
rho = 1.225; % Density of Air in kg/m3
u = 5; % air speed in m/s
Ambient_Temp = 30;
Water_Temp = 40;
RH = 40;
Water_FlowRate = 6; % in LPH
Water_FlowRate = Water_FlowRate/3600; % Conversion to kg/s
%Calculating Absolute Humidity for w_g()
abs = (6.112 * exp((17.67 * Ambient_Temp)/(Ambient_Temp + 243.5)) * RH * 2.1674)/((Ambient_Temp+273.15)) * (1/1000);
ma = rho * u * W * d; % mass flow rate of the air (per pan) in kg/s
Kel = 273.15; % Celsius to Kelvin conversion
patm = 101325; % atmospheric pressure in Pa
delX = 0.01; % Size of each step
max_iters = 1/delX; % Number of Iterations
cpv = 1864; % cp of vapor in J/kg-K
cpa = 1005; % cp of dry air in J/kg-K
cpl = 4182; % cp of water in J/kg-K
ifgw = 2500 * 1000; % Latent heat of water at 0C in J/kg
% Constants in saturation pressure definition
a1 = -5800.2206;
a2 = 1.3914993;
a3 = -0.048640239;
a4 = 0.41764768 * 10^-4;
a5 = -0.14452093 * 10^-7;
a6 = 6.5459673;
%% Arrays
w_g = 1:max_iters;
T_g = 1:max_iters;
T_l = 1:max_iters;
ml = 1:max_iters;
ig = 1:max_iters;
psat = 1:max_iters;
iv = 1:max_iters;
w_int = 1:max_iters;
h = 1:max_iters;
Lef = 1:max_iters;
cpg = 1:max_iters;
k = 1:max_iters;
w_g(1) = abs/rho; % in relation to absolute humidity calculated above (w_g) is dimensionless so one must divide it by rho so the kg/kg cancels out
T_g(1) = Ambient_Temp + Kel; % Temperature of the Gas
T_l(1) = Water_Temp + Kel; % Temperature of the liquid
ml(1) = Water_FlowRate; % Mass flow rate of the liquid in Kg/s
ig(1) = (cpa)*(T_g(1) - Kel) + (w_g(1) * (ifgw + (cpv * (T_g(1) - Kel)))); % Enthalpy of the moist air in Joules/Kg
%% Dependent terms
va = ((0.009392 * T_l(1)) - 1.237) * 10^-5; % dynamic viscosity




%
Dh = (2 * W * d)/(W + d); % Hydraulic Diameter 
Re = (u * Dh)/ va; % Reynolds number calculation




Pr = (-0.0002697 * T_l(1)) + 0.81; %Prandtl number




n = 0.4;
lambda = 0.027; % Conductivity of air in W/m-K
psat(1) = exp((a1/T_l(1)) + (a2) + (a3 * T_l(1)) + (a4 * (T_l(1))^2) + (a5 * (T_l(1))^3) + (a6 * log(T_l(1)))); % Saturation pressure of the water

iv(1) = ifgw + cpv*(T_l(1) - Kel); % enthalpy of the water vapour

w_int(1) = (0.622 * psat(1))/(patm-psat(1)); % Moisture content at the interface

if (T_l(1) < T_g(1))
    n = 0.3;
end
h(1) = 0.023*(Re^0.8)*(Pr^n)*(lambda/Dh); % Defining the constant h

Lef(1) = ((0.865)^(2/3))*(((w_int(1) + 0.622)/(w_g(1) + 0.622) - 1)/(log((w_int(1) + 0.622)/(w_g(1) + 0.622)))); % Defining the Lewis constant

cpg(1) = cpa + w_g(1) * cpv; % Defining the specific heat of moist air 

k(1) = h(1)/(cpg(1) * Lef(1)); % Definition of constant K
%% Runge-Kutta 4th order Method 
for ctr = 2:max_iters
    
    
    % Defining first terms
    c = ctr - 1;
    
    ki1 = k(c);
    wint1 = w_int(c);
    wg1 = w_g(c);
    Tl1 = T_l(c);
    Tg1 = T_g(c);
    iv1 = iv(c);
    h1 = h(c);
    ml1 = ml(c);
    ig1 = ig(c);
    
    pd = wint1 - wg1;
    td = Tl1 - Tg1;
    
    j1 = delX * (-1 * ki1 * A * pd);
    k1 = delX * (ki1 * A * pd / ma);
    l1 = delX * ((h1 * A * td / ma) + (iv1 * ki1 * A * pd / ma));
    o1 = delX * ((A / (ml1 * cpl)) * ((ki1 * pd * (cpl*(Tl1 - Kel) - iv1)) - (h1 * td)));
    
    % Defining second terms
    
    ml2 = ml1 + (j1/2);
    wg2 = wg1 + (k1/2);
    ig2 = ig1 + (l1/2);
    Tl2 = Tl1 + (o1/2);
    
    va = ((0.009392 * Tl2) - 1.237) * 10^-5; % dynamic viscosity
    
    Re = (u * Dh)/ va; % Reynolds number calculation
    Pr = (-0.0002697 * Tl2) + 0.81; %Prandtl number
    
    psat2 = exp((a1/Tl2) + (a2) + (a3 * Tl2) + (a4 * (Tl2)^2) + (a5 * (Tl2)^3) + (a6 * log(Tl2)));
    iv2 = ifgw + cpv*(Tl2 - Kel);
    
    wint2 = (0.622 * psat2)/(patm-psat2);
    Tg2 = ((ig2 - (wg2*ifgw))/(cpa + wg2*cpv)) + Kel;
    
    if (Tl2 < Tg2)
        n = 0.3;
    end
    
    h2 = 0.023*(Re^0.8)*(Pr^n)*(lambda/Dh);
    
    Lef2 = ((0.865)^(2/3))*(((wint2 + 0.622)/(wg2 + 0.622) - 1)/(log((wint2 + 0.622)/(wg2 + 0.622))));
    cpg2 = cpa + wg2 * cpv;
    
    ki2 = h2/(cpg2 * Lef2);
    
    pd2 = wint2-wg2; 
    td2 = Tl2-Tg2;
    
    j2 = delX * (-1 * ki2 * A * pd2);
    k2 = delX * (ki2 * A * pd2 / ma);
    l2 = delX * ((h2 * A * td2 / ma) + (iv2 * ki2 * A * pd2 / ma));
    o2 = delX * ((A / ml2 / cpl) * (ki2 * pd2 * (cpl*(Tl2-Kel) - iv2) - (h2 * td2)));
    
    % Defining third terms
    
    ml3 = ml1 + (j2/2);
    wg3 = wg1 + (k2/2);
    ig3 = ig1 + (l2/2);
    Tl3 = Tl1 + (o2/2);
    
    va = ((0.009392 * Tl3) - 1.237) * 10^-5; % dynamic viscosity
    Re = (u * Dh)/ va; % Reynolds number calculation
    Pr = (-0.0002697 * Tl3) + 0.81; %Prandtl number
    
    psat3 = exp((a1/Tl3) + (a2) + (a3 * Tl3) + (a4 * (Tl3)^2) + (a5 * (Tl3)^3) + (a6 * log(Tl3)));
    iv3 = ifgw + cpv*(Tl3 - Kel);
    
    wint3 = (0.622 * psat3)/(patm-psat3);
    Tg3 = ((ig3 - (wg3*ifgw))/(cpa + wg3*cpv)) + Kel;
    
    if (Tl3 < Tg3)
        n = 0.3;
    end
    
    h3 = 0.023*(Re^0.8)*(Pr^n)*(lambda/Dh);
    
    Lef3 = ((0.865)^(2/3))*(((wint3 + 0.622)/(wg3 + 0.622) - 1)/(log((wint3 + 0.622)/(wg3 + 0.622))));
    cpg3 = cpa + wg3 * cpv;
    
    ki3 = h3/(cpg3 * Lef3);
    
    pd3 = wint3-wg3; 
    td3 = Tl3-Tg3;
    
    j3 = delX * (-1 * ki3 * A * pd3);
    k3 = delX * (ki3 * A * pd3 / ma);
    l3 = delX * ((h3 * A * td3 / ma) + (iv3 * ki3 * A * pd3 / ma));
    o3 = delX * ((A / ml3 / cpl) * (ki3 * pd3 * (cpl*(Tl3-Kel) - iv3) - (h3 * td3)));
    
    % Defining 4th terms
    
    ml4 = ml1 + (j3/2);
    wg4 = wg1 + (k3/2);
    ig4 = ig1 + (l3/2);
    Tl4 = Tl1 + (o3/2);
    
    va = ((0.009392 * Tl4) - 1.237) * 10^-5; % dynamic viscosity
    Re = (u * Dh)/ va; % Reynolds number calculation
    Pr = (-0.0002697 * Tl4) + 0.81; %Prandtl number
    
    psat4 = exp((a1/Tl4) + (a2) + (a3 * Tl4) + (a4 * (Tl4)^2) + (a5 * (Tl4)^3) + (a6 * log(Tl4)));
    iv4 = ifgw + cpv*(Tl4 - Kel);
    
    wint4 = (0.622 * psat4)/(patm-psat4);
    Tg4 = ((ig4 - (wg4*ifgw))/(cpa + wg4*cpv)) + Kel;
    
    if (Tl4 < Tg4)
        n = 0.3;
    end
    
    h4 = 0.023*(Re^0.8)*(Pr^n)*(lambda/Dh);
    
    Lef4 = ((0.865)^(2/3))*(((wint4 + 0.622)/(wg4 + 0.622) - 1)/(log((wint4 + 0.622)/(wg4 + 0.622))));
    cpg4 = cpa + wg4 * cpv;
    
    ki4 = h4/(cpg4 * Lef4);
    
    pd4 = wint4-wg4; 
    td4 = Tl4-Tg4;
    
    j4 = delX * (-1 * ki4 * A * pd4);
    k4 = delX * (ki4 * A * pd4 / ma);
    l4 = delX * ((h4 * A * td4 / ma) + (iv4 * ki4 * A * pd4 / ma));
    o4 = delX * ((A / ml4 / cpl) * (ki4 * pd4 * (cpl*(Tl4-Kel) - iv4) - (h4 * td)));
    
    % Finding the value at i+1 
    
    ml(ctr) = ml(c) + ((j1 + 2 * j2 + 2 * j3 + j4)/6);
    w_g(ctr) = w_g(c) + ((k1 + 2 * k2 + 2 * k3 + k4)/6);
    ig(ctr) = ig(c) + ((l1 + 2 * l2 + 2 * l3 + l4)/6);
    T_l(ctr) = T_l(c) + ((o1 + 2 * o2 + 2 * o3 + o4)/6);
    
    if (Tl1 <= 0) || (Tl2 <= 0) || (Tl3 <= 0) || (Tl4 <= 0)
        
        disp("ERROR!");
        disp("1 : " + wint1 + " psat : " + psat(c) + " Tl : " + Tl1 + " o value : " + o1 + " TL : " + Tl1);
        disp("2 : " + wint2 + " psat : " + psat2 + " Tl : " + Tl2 + " o value : " + o2 + " TL : " + Tl2);
        disp("3 : " + wint3 + " psat : " + psat3 + " Tl : " + Tl3 + " o value : " + o3 + " TL : " + Tl3);
        disp("4 : " + wint4 + " psat : " + psat4 + " Tl : " + Tl4 + " o value : " + o4 + " TL : " + Tl4);
        
        disp(ctr);
        break;
    end
    
   % Updating dependent values 
   
    T_g(ctr) = ((ig(ctr) - (w_g(ctr)*ifgw))/(cpa + w_g(ctr)*cpv)) + Kel;
    
    va = ((0.009392 * T_l(ctr)) - 1.237) * 10^-5; % dynamic viscosity
    Re = (u * Dh)/ va; % Reynolds number calculation
    Pr = (-0.0002697 * T_l(ctr)) + 0.81; %Prandtl number
    
    psat(ctr) = exp((a1/T_l(ctr)) + (a2) + (a3 * T_l(ctr)) + (a4 * (T_l(ctr))^2) + (a5 * (T_l(ctr))^3) + (a6 * log(T_l(ctr)))); % Saturation pressure of the water
    iv(ctr) = ifgw + cpv*(T_l(ctr) - Kel); % enthalpy of the water vapour
    w_int(ctr) = (0.622 * psat(ctr))/(patm-psat(ctr)); % Moisture content at the interface
    if (T_l(ctr) < T_g(ctr))
    n = 0.3;
    end
    h(ctr) = 0.023*(Re^0.8)*(Pr^n)*(lambda/Dh); % Defining the constant h
    Lef(ctr) = ((0.865)^(2/3))*(((w_int(ctr) + 0.622)/(w_g(ctr) + 0.622) - 1)/(log((w_int(ctr) + 0.622)/(w_g(ctr) + 0.622)))); % Defining the Lewis constant
    cpg(ctr) = cpa + w_g(ctr) * cpv; % Defining the specific heat of moist air 
    k(ctr) = h(ctr)/(cpg(ctr) * Lef(ctr)); % Definition of constant K   
   
end
%%
L_array = L/max_iters:L/max_iters:L;
t_w = 1/10000;
t = L/ma;
diff = ml(1) - ml(end);
evap = diff * 3600;
dw = w_int - w_g;
O_AH = w_g(end)* 1.225;
O_Tg = T_g(end) - Kel;
O_RH = (1000 * O_AH * (273.15 + O_Tg)) / (6.112 * exp((17.67 * O_Tg)/(O_Tg + 243.5)) * 2.1674);
disp("Output RH is : " + O_RH + "%");

Best Answer

Related Solutions

MATLAB: How to write a code for iteration

MATLAB: My fluid simulation exploding to very large numbers when I take small step values

Related Question