Izračun enačbe

[matson]

V Excelu sem napisal serijo enačb, kjer vstavljam vrednosti in se vse samodejno računa. Sedaj pa sem pri tem postopku prišel do naslednje enačbe, ki je ne znam izračunati:

https://latex.codecogs.com/gif.latex?%5 ... ight%20%29


Poznam vse podatke, izračunati pa bi rabil Vr.

Kako bi se tega lotil v Excelu? Imam sicer tudi Matlab, a ga ne znam prav veliko.

P.S.: Kakor vidim, ni možno pripeti Excelove datoteke?


Bi mi lahko kdo popravil spodnjo kodo, da bo delovala na forumu?

\[\frac{P_{r}V_{r}}{T_{r}}=1+\frac{B}{V_{r}}+\frac{C}{V_{r}^{2}}+\frac{D}{V_{r}^{5}}+\frac{c_{4}}{T_{r}^{3}V_{r}^{2}}\left ( \beta +\frac{\gamma }{V_{r}^{2}} \right )exp\left ( -\frac{\gamma }{V_{r}^{2}} \right )\]

[shrink]

Ker gre za transcendentno enačbo, je nemogoče direktno (analitično) izraziti \(V_r\). Lahko pa se računa numerično, kar z Excelom ni možno, z Matlabom pa je.

[matson]

V Matlabu sem vstavil kodo:

syms Vr_R
eqn = Pr_0*Vr_R/Tr_0 == 1 + B_R/Vr_R + C_R/Vr_R^2 + D_R/Vr_R^5 + c4_R/(Tr_0^3*Vr_R^2)*(BETA_R + (GAMA_R/Vr_R^2)*exp(-GAMA_R/Vr_R^2)) ;
Vr_R=double(solve(eqn))

Rezultat je sicer pravi, vendar mi program javi:

Warning: Cannot solve symbolically. Returning a numeric approximation instead.

Tudi računa kakšno sekundo dlje, vendar me to ne moti.

Poizkušal sem s funkcijo fsolve, a ne znam.

[shrink]

No, očitno ti je Matlab, kljub poskusu simboličnega računanja, vrnil numeričen rezultat, kar je pač to, kar si iskal. Lahko direktno numerično računaš (s kakšno numerično funkcijo), a bo rezultat enak. Če pa to nujno potrebuješ, povej, pa ti bomo skušali pomagati s kodo.

[matson]

Gre za vmesni rezultat, a tisti končni je povsem pravi. Kakor sem gledal, se to da narediti z neko funkcijo fzero, ampak je kar malo zakomplicirano.

Ali je možno, da s to sedanjo kodo dobim napačen rezultat? Ali obstajata le dve varianti - pravi rezultat ali pa program sploh ne izračuna ničesar?

Bi pa rad poskusil z fzero (ali katero drugo), saj bo v nadaljnjem delu še katera podobna enačba, ki jo "solve" morda ne bo rešil....

[matson]

Sedaj sem malo poskušal z različnimi vrednostmi in mi pri nekaterih ne deluje:

Error using ^
Inputs must be a scalar and a square matrix.
To compute elementwise POWER, use POWER (.^) instead.

Error in KOMPLET_POSKUS (line 105)
eqn = Zr_R == 1 + B_R/Vr_R + C_R/Vr_R^2 + D_R/Vr_R^5 +
c4_R/(Tr_0^3*Vr_R^2)*(BETA_R + (GAMA_R/Vr_R^2)*exp(-GAMA_R/Vr_R^2)) ;

Zanima pa me, kako bi nastavil naslednje:

V Matlabu sem izračunal enačbo za 2 podani vrednosti (input =...). Sedaj pa bi rad, da mi ti dve vrednosti pobira avtomatsko in izračuna končne podatke.

Se pravi:
Pr, kjer so vrednosti od 0.01; 0.05; 0.1; 0.2; 0.4; 0.6; 0,8; 1; 1.2; 1.5; 2; 3; 5; 7; 10
Tr, kjer so vrednosti od 0.30; 0.35; 0.40 ...... 4.00 (korak med njimi pa ni vseskozi enak)

Izračunati pa bi mi moral vse možne kombinacije. Zaporedje pri tem ni pomembno, ker bom rezultate vstavil v tabelo. Tako, da bom lahko na koncu pogledal - pri Tr = 0.4 in Pr = 1, je rezultat enačbe "tak".

Bi kdo to znal nastaviti?

[shrink]

Navedi še vrednosti za koeficiente, da preizkusim kodo v Matlabu.

P.S. Navedi tudi celoten nabor vrednosti za \(T_r\).

[matson]

Sem prilepil kar kodo (prvi del je za osnovni fluid, drugi pa za referenčni).

Velja, da je Tr_0 = T_0/Tc_0.
Velja, da je Tr_R = T_R/Tc_R.
Velja, da je Pr_0 = P_0/Pc_0.
Velja, da je Pr_R = P_R/Pc_R.

Sedaj je nastavljeno, da vstavljam T_0 in P_0 (vrstici 23 in 24) program pa mi računa Tr_0 in Pr_0 (vrst. 31 in 35). Enako potem sledi v drugem delu kode za referenčni fluid. Če bom vstavljal različne Pr_0 in Tr_0 (ter Pr_R in Tr_R), potem teh vrstic ne bom rabil (vrstice 22 - 35 in 70 - 83).

Vrednosti za Tr so: 0.30; 0.35; 0.40; 0.45; 0.50; 0.55; 0.60; 0.70; 0.75; 0.80; 0.85; 0.90; 0.93; 0.95; 0.97; 0.98; 0.99; 1.00; 1.01; 1.02; 1.05; 1.1; 1.15; 1.20; 1.30; 1.40; 1.50; 1.60; 1.70; 1.80; 1.90; 2.00; 2.20; 2.40; 2.60; 2.80; 3.00; 3.50; 4.00
Vrednosti za Pr so: 0.010; 0.050; 0.100; 0.200; 0.400; 0.600; 0.800; 1.000; 1.200; 1.500; 2.000; 3.000; 5.000; 7.000; 10.000

Rabim pa vse kombinacije, da bom imel na koncu tabelo (matriko). To bom kopiral v Excel in naredil tabelo.

%% PROGRAM ZA IZRAČUN FAKTORJA STISLJIVOSTI PO METODI LEE-KESLER
%%

'_0 - OSNOVNI FLUID'
'_R - REFERENČNI FLUID'

%%
% *'OSNOVNI FLUID'*

%Konstante za osnovni fluid
b1_0=0.1181193; b2_0=0.265728; b3_0=0.154790; b4_0=0.03023; c1_0=0.0236744; c2_0=0.0186984; c3_0=0; c4_0=0.042724; d1_0=0.155488*10^-4;
d2_0=0.623689*10^-4; BETA_0=0.65392; GAMA_0=0.060167;

%IZBERI ŽELJENI OSNOVNI FLUID
%Kritične vrednosti za osnovni fluid - METAN
Tc_0=190.564; Pc_0=4599000; Vc_0=0.0986; OMEGAc_0=0.286;
%Kritične vrednosti za osnovni fluid - VODA
%Tc_0=647.096; Pc_0=22064000; Vc_0=0.0559472; OMEGAc_0=0.229;
%Kritične vrednosti za osnovni fluid - ETAN
%Tc_0=305.32; Pc_0=4872000; Vc_0=0.1455; OMEGAc_0=0.09493;

%Podatki za osnovni fluid, T[K], p[Pa]
T_0=input('VSTAVI VREDNOST TEMPERATURE OSNOVNEGA FLUIDA V [K]:')
P_0=input('VSTAVI VREDNOST TLAKA OSNOVNEGA FLUIDA V [Pa]:')
%T_0=288
%P_0=101325

%Izračun reduciranih vrednosti za osnovni fluid
syms Tr_0
eqn = Tr_0 == T_0/Tc_0 ;
Tr_0=double(solve(eqn))

syms Pr_0
eqn = Pr_0 == P_0/Pc_0 ;
Pr_0=double(solve(eqn))

%Izračun konstant B_0, C_0, D_0 (osnovni fluid)
syms B_0
eqn = B_0 == b1_0 - b2_0/Tr_0 - b3_0/Tr_0^2 - b4_0/Tr_0^3 ;
B_0=double(solve(eqn))

syms C_0
eqn = C_0 == c1_0 - c2_0/Tr_0 + c3_0/Tr_0^3 ;
C_0=double(solve(eqn))

syms D_0
eqn = D_0 == d1_0 + d2_0/Tr_0 ;
D_0=double(solve(eqn))

%Izračun Vr_0 za osnovni fluid
syms Vr_0
eqn = Pr_0*Vr_0/Tr_0 == 1 + B_0/Vr_0 + C_0/Vr_0^2 + D_0/Vr_0^5 + c4_0/(Tr_0^3*Vr_0^2)*(BETA_0 + (GAMA_0/Vr_0^2)*exp(-GAMA_0/Vr_0^2)) ;
Vr_0=double(solve(eqn))

%Izračun Z0_0
syms Z0_0
eqn = Z0_0 == 1 + B_0/Vr_0 + C_0/Vr_0^2 + D_0/Vr_0^5 + c4_0/(Tr_0^3*Vr_0^2)*(BETA_0 + (GAMA_0/Vr_0^2)*exp(-GAMA_0/Vr_0^2)) ;
Z0_0=double(solve(eqn))


%%
% *'REFERENČNI FLUID - N-OKTAN'*

%Konstante za referenčni fluid
b1_R=0.2026579; b2_R=0.331511; b3_R=0.027655; b4_R=0.203488; c1_R=0.0313385; c2_R=0.0503618; c3_R=0.016901; c4_R=0.041577; d1_R=0.48736*10^-4;
d2_R=0.0740336*10^-4; BETA_R=1.226; GAMA_R=0.03754
%Kritične vrednosti za referenčni fluid
Tc_R=562,05; Pc_R=4895000; Vc_R=0.256; OMEGAc_R=0.3978

%Podatki za referenčni fluid (n-oktane), T[K], p[Pa]
T_R=input('VSTAVI VREDNOST TEMPERATURE REFERENČNEGA FLUIDA V [K]:')
P_R=input('VSTAVI VREDNOST TLAKA REFERENČNEGA FLUIDA V [Pa]:')
%T_R=288
%P_R=101325

%Izračun reduciranih vrednosti za referenčni fluid
syms Tr_R
eqn = Tr_R == T_R/Tc_R ;
Tr_R=double(solve(eqn))

syms Pr_R
eqn = Pr_R == P_R/Pc_R ;
Pr_R=double(solve(eqn))

%Izračun konstant B_R, C_R, D_R (referenčni fluid)
syms B_R
eqn = B_R == b1_R - b2_R/Tr_R - b3_R/Tr_R^2 - b4_R/Tr_R^3 ;
B_R=double(solve(eqn))

syms C_R
eqn = C_R == c1_R - c2_R/Tr_R + c3_R/Tr_R^3 ;
C_R=double(solve(eqn))

syms D_R
eqn = D_R == d1_R + d2_R/Tr_R ;
D_R=double(solve(eqn))

%Izračun Vr za referenčni fluid - ali je Vr od osnovnega ali od referenčnega?????? Pomoje referenčni!!!
syms Vr_R
eqn = Pr_0*Vr_R/Tr_0 == 1 + B_R/Vr_R + C_R/Vr_R^2 + D_R/Vr_R^5 + c4_R/(Tr_0^3*Vr_R^2)*(BETA_R + (GAMA_R/Vr_R^2)*exp(-GAMA_R/Vr_R^2)) ;
Vr_R=double(solve(eqn))

%Izračun Zr_R
syms Zr_R
eqn = Zr_R == 1 + B_R/Vr_R + C_R/Vr_R^2 + D_R/Vr_R^5 + c4_R/(Tr_0^3*Vr_R^2)*(BETA_R + (GAMA_R/Vr_R^2)*exp(-GAMA_R/Vr_R^2)) ;
Zr_R=double(solve(eqn))


%% 'Izračun Z - faktor stisljivosti'*
syms Z
eqn = Z == Z0_0 + OMEGAc_0/OMEGAc_R * (Zr_R - Z0_0) ;
Z=double(solve(eqn))
format long

[shrink]

Tr_0 in Pr_0 lahko definiraš kot polji vrednosti, npr. tako:

Koda: Izberi vse

Tr_0=[0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.70, 0.75, 0.80, 0.85, 0.90, 0.93, 0.95, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.05, 1.1, 1.15, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00, 2.20, 2.40, 2.60, 2.80, 3.00, 3.50, 4.00]
Pr_0=[0.010, 0.050, 0.100, 0.200, 0.400, 0.600, 0.800, 1.000, 1.200, 1.500, 2.000, 3.000, 5.000, 7.000, 10.000]

za avtomatizirano reševanje za različne vrednosti pa lahko uporabiš dvojno 'for' zanko, npr. tako:

Koda: Izberi vse

for j =1:39
    for i = 1:15
        syms Vr_R
        solve(-Pr_0(i)*Vr_R/Tr_0(j) + 1 + B_R/Vr_R + C_R/Vr_R^2 + D_R/Vr_R^5 + c4_R/(Tr_0(j)^3*Vr_R^2)*(BETA_R + 	
        GAMA_R/Vr_R^2)*exp(-GAMA_R/Vr_R^2), Vr_R)
    end
end

pri čemer je Tr_0 nadomeščen s Tr_0(j) in Pr_0 s Pr_0(i).

Reševanje na gornji način s funkcijo 'solve' daje obilo kompleksnih in le za vzorec realnih rešitev. Uporabil sem obliko enačbe, ki si jo produciral z LaTeX kodo na vrhu (tvoj zapis v Matlabovi kodi je rahlo drugačen: oklepaji namreč niso na pravih mestih). Preveri še enkrat, če si na vrhu zapisal pravilno enačbo oz. zakaj ni v skladu s tvojim zapisom. Morda je to vzrok za kompleksne rešitve, morda pa tudi algoritem nepravilno konvergira k rešitvam.

P.S. Sam uporabljam eno starejšo različico Matlaba, ki ne podpira sintakse za funkcijo 'solve', kot jo sam uporabljaš. Za numerično reševanje v novih različicah sicer obstaja funkcija 'vpasolve':

http://www.mathworks.com/help/symbolic/vpasolve.html

[shrink]

P.S.2 Sem preveril še Maple-u: numerične rešitve so popolnoma enake.

[matson]

Ja, tista v Latexu je prava. Na koncu sem napačno postavil zaklepaj. Še enkrat sem preveril za določen Tr in Pr. Rezultat je enak. Očitno je majhna vrednost.

Super. Tole sem vstavil v mojo kodo in deluje. Praviš, da si preveril rezultate. Sedaj mi računa kodo že kar nekaj časa.

Rad pa bi naredil, da mi vrže ven matriko oz. tabelo, saj bom na koncu moral iz teh rezultatov oblikovati tabelo (najlažje v Excelu). Tudi kopiranje matrike, bo najbrž imelo svoj čar.

Koda: Izberi vse

a = [1, 3, 6, 7, 10, 12]
b = [2, 3, 8, 17, 25, 44]

X = a + b
X = a .* b
X = bsxfun(@plus, a', b)

Sedaj sem tej tvoji kodi na koncu dodal še:

Koda: Izberi vse

X = bsxfun(@plus, Tr_0', Pr_0)

Sem pa v Matlabu od prejšnjega tedna, tako da vse delam povsem na pamet. :/

To niti ni tako nujno, a vseeno - Ali se da narediti tako, da mi ne računa vseh rezultatov posamezno, temveč vrže ven matriko?

Mi pa pri vsakem rezultatu še vedno napiše:

Koda: Izberi vse

Warning: Cannot solve symbolically. Returning a numeric approximation
instead. 

Potem moram vse te rezultate vnesti v enačbo za Z0_0 in te prikazati v tabeli. Enako potem z Zr_0. Tabelo bolj rabim za te končne Z-je.

[matson]

Spodnja koda deluje in mi na koncu izriše matriko rezultatov (nisem jih preveril).

Koda: Izberi vse

Tr_0=[1.34, 1.515789]
Pr_0=[0.0220319, 1.345, 1.5]

for j =1:2
    for i = 1:3
        syms Vr_0
        solve(-Pr_0(i)*Vr_0/Tr_0(j) + 1 + B_0/Vr_0 + C_0/Vr_0^2 + D_0/Vr_0^5 + c4_0/(Tr_0(j)^3*Vr_0^2)*(BETA_0 + GAMA_0/Vr_0^2)*exp(-GAMA_0/Vr_0^2), Vr_0)
    end
end
Vr_0 = bsxfun(@plus, Tr_0', Pr_0)

Sedaj pa bi moral to povezati s kodo za Z0_0 in izrisati tabelo (matriko) rezultatov. Če dodam spodnjo kodo, mi ne deluje:

Koda: Izberi vse

%Izračun Z0_0
for j =1:6
    syms Z0_0
    solve(-Z0_0 + 1 + B_0/Vr_0(j) + C_0/Vr_0(j)^2 + D_0/Vr_0(j)^5 + c4_0/(Tr_0^3*Vr_0(j)^2)*(BETA_0 + (GAMA_0/Vr_0(j)^2))*exp(-GAMA_0/Vr_0(j)^2), Z0_0)
end
Z0_0 = bsxfun(@plus, Tr_0', Pr_0)

Error

Koda: Izberi vse

Error using  ^ 
Inputs must be a scalar and a square matrix.
To compute elementwise POWER, use POWER (.^) instead.

Error in dfgherthe (line 30)
    solve(-Z0_0 + 1 + B_0/Vr_0(j) + C_0/Vr_0(j)^2 + D_0/Vr_0(j)^5 +
    c4_0/(Tr_0^3*Vr_0(j)^2)*(BETA_0 +
    (GAMA_0/Vr_0(j)^2))*exp(-GAMA_0/Vr_0(j)^2), Z0_0)

Poskusil sem dodati še eno for zanko, a tudi ne gre.

[shrink]

Ne deluje ti zato, ker je Vr_0 matrika, torej v splošnem dvodimenzionalno polje Vr_0(i,j), ti pa si ga obravnaval kot enodimenzionalno polje Vr_0(j). Moraš torej klicati Vr_0(i,j), seveda spet v dvojni 'for' zanki.

[matson]

Sedaj sem opazil. V naslednjo enačbo dejansko ne vnašam Pr_0 in Tr_0, temveč le Tr_0 in matriko Vr_0. Kako pa naredim, da kličem to matriko in obenem še Tr_0?

[shrink]

Tako kot prej: s Tr_0(j).