How to improve accuracy issues with solve function or just code itself?

Hi, I'm writing a code to find the viscosity of a fluid given a user input temperature (T) and pressure (P) (where the pressure has to be between 0.01 and 100). This is my entire script:

clear
clc

%formula 1
A = [0.128517, -1.00, -0.709661, 0.662534, -0.197846, 0.00770147];
B = [0.465601, 1.26469, -0.511425, 0.2746];
A05 = 2.60661;
pstar = 314.4; 
Tstar = 132.5; 
H = 6.1609*10^-6; 

%formula 2
MWair = 0.029; % 0.79 oxygen with 0.21 nitrogen
Bstar = 20.5*10^-6;
Tb = 340;
b = [-0.022243, 2.16059, -1.94783, -0.1704, -0.0223299, -0.00168785];
c = [1.19665, 4.29631, -5.52109, 2.50287];
R = 8.314;

T=input('Input Temperature (K): ');
P=input('Input Pressure (between 0.01 and 1000 MPa): ');

Tr = T/Tstar;

while (P<0.01 || P>1000)
    disp(' ');
    P = input('Invalid pressure. Please input pressure between 0.01 and 1000 MPa: ');
end

%calculating initial viscosity
n0=A05*Tr^(0.5);
for i=1:6
   n0=n0+A(i)*Tr^(2-i);
end

%density calculation
tau = T/Tb;
Btau = 0;
Ctau = 0;
for i=1:6
    Btau = Btau + b(i)*tau^(2-i);
end
for j=1:4
    Ctau = Ctau + c(j)*tau^(1-j);
end

Btau = Btau*Bstar;
Ctau = Ctau*Bstar^2;

syms pbar;
pb = solve(1 + Btau*pbar + Ctau*pbar^2 - P/(pbar*R*T), pbar,'Real',true);
pmol = vpa(pb);

p = pmol*MWair;
pr = p/pstar;
deltan = 0;

for i=1:4
    deltan = deltan + B(i)*pr^i;
end

n = H*(n0 + deltan);

display('Viscosity given pressure and temperature (Pa*s): ')
display(n)

It is required that I related Z = 1 + B(pbar) + C(pbar)^2 and Z = P/(pbar*R*T). It is also necessary that the calculated viscosity (which I defined as variable 'n') is accurate within three significant figures. I am having trouble reaching accuracy with high temperatures.

side note:For loops were made to achieve this general formula sum(AiXi^i) where Ai is any constant (A,B,b,c) at value i and Xi is an observed variable (pr and tau) at point i.

   T         P               n
  250     0.01          1.6040E-05 accurate
  300     0.5          1.8610E-05    |
  350     0.1          2.0900E-05 relatively accurate
  400      1            2.3170E-05 not so accurate  
  ..... skipped a few .....
  2000     100          7.0690E-05 horribly inaccurate

What can I do to improve the accuracy of my values? Are my accuracy errors due to the solve function?

Answer 1

[Walter Roberson]

You could convert all of your constants such as in "A" to be sym():

A = [sym('0.128517'), sym('-1.00'), sym(-0.709661'), sym('0.662534'), sym('-0.197846'), sym('0.00770147')]

including

T = sym(input('Input Temperature (K): ', 's'));

however, leave P as numeric until you have finished testing to see if it is in range, and then sym() it immediately after.

Recode

n0=A05*Tr^(0.5);

to

n0=A05*sqrt(Tr);

At the point you have

pmol = vpa(pb);

just assign

pmol = pb;

so that pmol stays symbolic.

Leave everything symbolic until finally what you do is

format long g
disp(double(n))

If you need higher precision than what you then get as a result, you will need to use

char(vpa(n, NumberOfDigits))    %eg 25 for NumberOfDigits

And if that shows inaccuracies then the problem is either in the calculation formula, or in insufficient accuracy of your constants.

Remember, when you have

MWair = 0.029

then any formula involving MWair can have at most 2 digits of precision since MWair is only 2 digits of precision. Likewise for the -1.00 that appears in "A", only two digits of precision.

[Jon]

Thank you for your response! Everything you said makes perfect sense. Unfortunately, for this project, I found that the largest problem to my code was that I did not convert the units of my pressure to the correct order of magnitude.

Nevertheless, thank you for taking the time to read my lengthy problem. Curious though, what does the format long g and disp(double(n)) do to symbolic values? Does it convert all the symbolic answers into single numerical answers?