Matlab problem

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?

Jaz bi se tega lotil takole: kreiraš dva vektorja, prvi vsebuje vse a-je, drugi vse b-je.
Nato bi z FOR zanko šel skozi prvi vektor a[i] za vse i-je, vgnezdeno v to pa še FOR zanko ki gre skozi drugi vektor b[j] za vse j-je. Znotraj te druge zanke bi izvajal tvoj ukaz (X = a + b).

Če uporabiš indeksa s katerimi se drajsaš skozi oba zgornja vektorja (recimo jima i in j)tudi kot definicijo lokacije v matriki za izhodne podatke X[i,j], boš dobil matriko z rezultati.
Srečno.

če sta a in b vektroja ne potrebuje for zank.
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) % za poljubne operatorje
...

V Octave moj primer deluje...

Vzami primer of bluefish in spremeni

X = i + j;
v
X = i + j + T + 20;

vendar se v matlabu probaj čimbolj izogibati zankam (for, while)...
zato lahko uporabiš

X = p(:) + q(:) + T + 20

namesto dveh for zank

Recimo kasneje bom imel tole enačbo. Odebeljene vrednosti so izračunane že od prej. Podčrtane vstavljam. Končni rezultat bo Vr.

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

Kot je rekel TEDY, zadnja vrstica je odveč (oz ni odveč, le ne prišteva še T-ja in fiksnega števila 20) Sicer je pa prvi X ki ga dobiš pravilen.

Končno enačbo premeči tako, da boš imel na levi samo še Vr_0 (Vr_0 == klobasa na desni * ( Tr_0/Pr_0) ), preveri mojo logiko ker sem iz vaje.

Poskusi kaj dobiš če v desni strani premetane enačbe tvoje podčrtane parametre zamenjaš z ustreznimi klici na vektorje Pr_0(:), Tr_(:). Če dobiš preveč rezultatov (ker izračunalo kombinacijo 3 spremenljivk - Tr_0 boš imel na desni strani dvakrat), se vrni na rešitev ki ti jo je podal Bluefish in tisti njegov X = i + j zamenjaj s svojo formulo.

Tole je, kar mi je uspelo do sedaj. V to enačbo sem vstavljal Tr_0 in Pr_0.

Tr_0=[1, 1.5]
Pr_0=[0.5, 1.3, 2]

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)

Končni rezultat je matrika:

Vr_0 =

1.362031900000000 2.685000000000000 2.840000000000000
1.537820900000000 2.860789000000000 3.015789000000000

Sedaj pa moram v naslednjo enačbo vstaviti Vr_0 (rezultati matrike) in spet Tr_0 iz prvega dela. Dobiti moram neko matriko, da bo lepo pregledno.

%Izračun Z0_0
for j = Vr_0
for i = 1:3
syms Z0_0
solve(-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)
end
end
Z0_0 = bsxfun(@plus, Tr_0', Pr_0)

Na tak način, kot zgoraj... Kako sploh pokličem matriko?

Res? Ker jih v Workspace vseskozi vidim. Tudi izračune za Vr_0, ki sem jih omenil zgoraj. V Workspace-u 2x kliknem na matriko in odpre mi tabelo rezultatov (tako kot v Excelu).

Kako pa bi lahko tole transcendentno enačbo rešil drugače?

Oz. kako rešim drugi del, kjer vstavljam - kličem matriko Vr_0 in vrednosti Tr_0 (tiste od 0.3 pa do 4.00)?

Tale koda deluje odlično. Teh funkcij ne poznam in ne vem, kaj počnem. Sem se pa malo zafrknil, ker sem prej računal te kontante B, C in D le za eno spremenljivko. Sedaj bi moral v kodo vstaviti enačbe vseh treh, ki pa nabira vrednosti le Tr_0.

Verjetno se to da narediti s temi zankami? Vse tiste b1, b2 itd pa so standardi. Tele 3 pa je treba izračunati vsakič posebej.

B = double(solve(-B_0 + b1_0 - b2_0/Tr_0(j) - b3_0/Tr_0(j)^2 - b4_0/Tr_0(j)^3))
C = double(solve(-C_0 + c1_0 - c2_0/Tr_0(j) + c3_0/Tr_0(j)^3))
D = double(solve(-D_0 + d1_0 + d2_0/Tr_0(j)))

Ali sem na pravi poti s spodnjo kodo?

Z0 = []; 
    for j = 1:length(Tr_0)
        B = double(solve(-B_0 + b1_0 - b2_0/Tr_0(j) - b3_0/Tr_0(j)^2 - b4_0/Tr_0(j)^3))
        C = double(solve(-C_0 + c1_0 - c2_0/Tr_0(j) + c3_0/Tr_0(j)^3))
        D = double(solve(-D_0 + d1_0 + d2_0/Tr_0(j)))
    for i = 1:length(Pr_0)
            Vr = double(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));
         
        %Izračun Z0_0
        Z0(end+1) = double(solve(-Z0_0 + 1 + B_0/Vr + C_0/Vr^2 + D_0/Vr^5 + c4_0/(Tr_0(j)^3*Vr^2)*(BETA_0 + (GAMA_0/Vr^2))*exp(-GAMA_0/Vr^2), Z0_0));
    end
end
 
disp(Z0)

Kaj pomeni (end+1)?

%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;

%clear all
Tr_0=[1.34, 1.515789]
Pr_0=[0.0220319, 1.345, 1.5]
syms B_0 C_0 D_0 Vr_0 Z0_0
 
Z0 = []; 
    for j = 1:length(Tr_0)
        B = (-B_0 + b1_0 - b2_0/Tr_0(j) - b3_0/Tr_0(j)^2 - b4_0/Tr_0(j)^3);
        C = (-C_0 + c1_0 - c2_0/Tr_0(j) + c3_0/Tr_0(j)^3);
        D = (-D_0 + d1_0 + d2_0/Tr_0(j));
    for i = 1:length(Pr_0)
            Vr = double(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));
          
        %Izračun Z0_0
        Z0(end+1) = double(solve(-Z0_0 + 1 + B_0/Vr + C_0/Vr^2 + D_0/Vr^5 + c4_0/(Tr_0(j)^3*Vr^2)*(BETA_0 + (GAMA_0/Vr^2))*exp(-GAMA_0/Vr^2), Z0_0));
    end
end
 disp(Z0)

Error:

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

Error in dfgherthe (line 51)
        Z0(end+1) = double(solve(-Z0_0 + 1 + B_0/Vr + C_0/Vr^2 + D_0/Vr^5
        + c4_0/(Tr_0(j)^3*Vr^2)*(BETA_0 +
        (GAMA_0/Vr^2))*exp(-GAMA_0/Vr^2), Z0_0));

PRINCIP NAJ BI BIL TAK:
1. Konstante (b1, b2...gama, beta) in Tr (od 0.3 pa do 4) vstavim v enačbe za B, C in D.
2. B, C, D + Tr (od 0.3 pa do 4) + Pr (od 0.05 do 10) vstavljam za enačbo Vr in dobim matriko rezultatov za Vr. Matrika, kjer so na y osi vstavljene vrednosti Tr, na x pa Pr.
3. Vse vstavim v zadnjo enačbo za Z0, kjer dobim matriko rezultatov Z0. Matrika, kjer so na y osi vstavljene vrednosti Tr, na x pa Pr. (to je edina matrika, ki jo dejansko rabim videti)

Je pa tole postalo zame malo prekomplicirano. Upam, da bo to na koncu vse zmešano.

Sem poskusil in deluje, vendar le za določene vrednosti Tr in Pr. Verjetno je fora v tisti transcendenti enačbi za Vr. V nekaterih delih tabele dobim pravilne rezultate (to velja za 85-90% tabele). Ko pa se pomikam proti manjšim vrednostim Tr (< 0,5) in Pr (< 0,4), pa so le še nekateri rezultati pravilni, na koncu pa so vsi povsem napačni.

Bolj pa me zanima naslednje. V nadaljevanju bom moral uporabiti to tabelo Z-jev za izračun različnih veličin. Tam bodo nastopali le Z-ji ter Tr in Pr (isti kot prej).

Npr.: Tole je enačba za izračun fugativnosti. Vsakič bom moral poračunati E (vseskozi enaka enačba) in nato vstaviti v veličino (5 jih imam, ena ima še nekaj zraven ostale so na isti način). To so potem končni rezultati. Vstavljati rabim dobljene vrednosti Vr in Z.

fugativnost = ln (f/p) = Z - 1 - ln(Z) + B/(2*Vr^2) + D/(5*Vr^5) + E

E = c4/(2*Tr^3*GAMA) (BETA + 1 - (BETA + 1 + GAMA/Vr^2)*exp(-GAMA/Vr^2))

Najbrž bi moral potem odpreti nov .m dokument, poklicati vrednosti Z-jev ali pa Vr-jev in potem napisati kodo za izračun npr. fugativnosti. Če ne bom znal rešiti tistega prvega problema, potem bi to lahko naredil za del, ki deluje avtomatično, ostalo pa spet poračunal kako drugače. Ali bom lahko iz te kode klical tudi Vr-je?

Druače je pa tole super koda!

Kako nastavim spodnjo kodo? Izračunane B, C in D vstavljam naprej v enačbi za Vr in Z0, kar je bilo prej v redu. Sedaj sem dodal še da se rezultati vstavijo v E in nato v F0, a mi jih ne bere.

%REDUCIRANA FUGATIVNOST ln(f/p)
syms Vr_0 Z0_0 E F0
F0 = []; 
    for j = 1:length(Tr_0)
        B_0 = (b1_0 - b2_0/Tr_0(j) - b3_0/Tr_0(j)^2 - b4_0/Tr_0(j)^3);
        C_0 = (c1_0 - c2_0/Tr_0(j) + c3_0/Tr_0(j)^3);
        D_0 = (d1_0 + d2_0/Tr_0(j));
    for i = 1:length(Pr_0)
            Vr = double(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));
            Z0 = double(solve(-Z0_0 + 1 + B_0/Vr + C_0/Vr^2 + D_0/Vr^5 + c4_0/(Tr_0(j)^3*Vr^2)*(BETA_0 + (GAMA_0/Vr^2))*exp(-GAMA_0/Vr^2), Z0_0));
            E = double(solve(-E + c4_0/(2*Tr_0^3) * (BETA + 1 - (BETA + 1 + GAMA/Vr^2)* exp(-GAMA/Vr^2))), E);
        %Izračun F0
        F0(end+1) = double(solve(-Z0_0 + 1 + B_0/Vr + C_0/Vr^2 + D_0/Vr^5 + c4_0/(Tr_0(j)^3*Vr^2)*(BETA_0 + (GAMA_0/Vr^2))*exp(-GAMA_0/Vr^2), F0));
    end 
end
disp(F0)

Res je, pa tudi tiste GAME in BETE niso bile ok.

Malo sem popravil, sedaj pa dobim nek čuden error (line 40 je v tem zapisu vrstica 18):

Error using sym/double
Too many input arguments.

Error in fugativnost (line 40)
            E = double(solve(-E + c4_0/(2*Tr_0(j)^3) * (BETA_0 + 1 -
            (BETA_0 + 1 + GAMA_0/Vr^2)* exp(-GAMA_0/Vr^2))), E);

Koda je taka:

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;

%clear all
Tr_0=[0.85, 0.90, 0.93, 0.95]
Pr_0=[0.800, 1.000, 1.200, 1.500]

%REDUCIRANA FUGATIVNOST ln(f/p)
syms Vr_0 Z0_0 E F0
F0 = []; 
    for j = 1:length(Tr_0)
        B_0 = (b1_0 - b2_0/Tr_0(j) - b3_0/Tr_0(j)^2 - b4_0/Tr_0(j)^3);
        C_0 = (c1_0 - c2_0/Tr_0(j) + c3_0/Tr_0(j)^3);
        D_0 = (d1_0 + d2_0/Tr_0(j));
    for i = 1:length(Pr_0)
            Vr = double(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));
            Z0 = double(solve(-Z0_0 + 1 + B_0/Vr + C_0/Vr^2 + D_0/Vr^5 + c4_0/(Tr_0(j)^3*Vr^2)*(BETA_0 + (GAMA_0/Vr^2))*exp(-GAMA_0/Vr^2), Z0_0));
            E = double(solve(-E + c4_0/(2*Tr_0(j)^3) * (BETA_0 + 1 - (BETA_0 + 1 + GAMA_0/Vr^2)* exp(-GAMA_0/Vr^2))), E);
        %Izračun F0
        F0(end+1) = double(solve(-Z0 + 1 + B_0/Vr + C_0/Vr^2 + D_0/Vr^5 + c4_0/(Tr_0(j)^3*Vr^2)*(BETA_0 + (GAMA_0/Vr^2))*exp(-GAMA_0/Vr^2), F0));
    end 
end
disp(F0)

Vidim, na koncu morata biti dva. Tole sem dal pa na isti način kot prej, a mi javi:

Error using solve>getEqns (line 422)
The input contains an empty equation or variable.

Error in solve (line 227)
[eqns,vars,options] = getEqns(varargin{:});

Error in fugativnost (line 42)
            F0 (end+1) = double(solve(-F0 + Z0 + 1 - log(Z0) + B_0/Vr +
            C_0/(2*Vr^2) + D_0/(5*Vr^5) + E, F0));

%REDUCIRANA FUGATIVNOST ln(f/p)
syms Vr_0 Z0_0 E F0
F0 = []; 
    for j = 1:length(Tr_0)
        B_0 = (b1_0 - b2_0/Tr_0(j) - b3_0/Tr_0(j)^2 - b4_0/Tr_0(j)^3);
        C_0 = (c1_0 - c2_0/Tr_0(j) + c3_0/Tr_0(j)^3);
        D_0 = (d1_0 + d2_0/Tr_0(j));
    for i = 1:length(Pr_0)
            Vr = double(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));
            Z0 = double(solve(-Z0_0 + 1 + B_0/Vr + C_0/Vr^2 + D_0/Vr^5 + c4_0/(Tr_0(j)^3*Vr^2)*(BETA_0 + (GAMA_0/Vr^2))*exp(-GAMA_0/Vr^2), Z0_0));
            E = double(solve(-E + c4_0/(2*Tr_0(j)^3) * (BETA_0 + 1 - (BETA_0 + 1 + GAMA_0/Vr^2)* exp(-GAMA_0/Vr^2)), E));
            %Izračun F0
            F0 (end+1) = double(solve(-F0 + Z0 + 1 - log(Z0) + B_0/Vr + C_0/(2*Vr^2) + D_0/(5*Vr^5) + E, F0));
    end
end
disp(F0)

Sedaj se s temle ubadam že nekaj časa. A je kaj zaradi tega (end+1)?

Je pa v Matlabu funkcija ln enaka kar log? Tam bi sicer moralo biti ln(Z0).

Ja, v Matlabu je log uporabljen za naravni logaritem, log10 pa za desetiškega.