Problem z dokładnością obliczeń

Witam wszystkich. Mój problem jest taki, że wyniki obliczeń z programu obarczone są sporym błędem. Mam porównanie do wyników z Mathcada, które są na pewno prawidłowe. Oto program:

#include <iostream>
#include <cmath>
#include <iomanip>
#include <vector>

using namespace std;

//prototypy funkcji
double ObliczCisnienieNasycenia(double &t_atm);
vector<double> ObliczWektorX1(double &p_atm, double &p_nas, 
	double &f, vector<double> ps, vector<double> a1);

int main()
{
	double p_atm, t_atm, f, p_nas, x_p;
	
	vector<double> ps;		//wektor współczynników pow. suchego
	ps.reserve(9);
	ps.push_back (1.009950160E+04);
	ps.push_back (-1.968275610E+02);
	ps.push_back (5.009155110E+00);
	ps.push_back (-5.761013730E-03);
	ps.push_back (1.066859930E-05);
	ps.push_back (-7.940297970E-09);
	ps.push_back (2.185231910E-12);
	ps.push_back (-1.767967310E+02);
	ps.push_back (-3.921504225E+00);

	vector<double> a1;		//wektor współczynników wody
	a1.reserve(9);
	a1.push_back (-3.947960830E+04);
	a1.push_back (5.755731020E+02);
	a1.push_back (9.317826530E-01);
	a1.push_back (7.222712860E-03);
	a1.push_back (-7.342557370E-06);
	a1.push_back (4.955043490E-09);
	a1.push_back (-1.336933246E-12);
	a1.push_back (-3.303974310E+04);
	a1.push_back (1.724205775E+01);

	t_atm = 288.15;
	f = 60;

	p_nas = ObliczCisnienieNasycenia(t_atm);
	cout << "Cisnienie nasycenia dla 15 *C wynosi: " << p_nas << " Pa" << endl;
	
	vector<double> X1;
	X1 = ObliczWektorX1(p_atm, p_nas, f, ps, a1);
	
	for(int i = 0; i < 9; i++) {
	cout <<  X1[i] << endl;
	}

	system("pause");
	return 0;
}

double ObliczCisnienieNasycenia(double &t_atm) {

	double V, A, B, C, p_sat;
	const double p = 1.0E06, t = 1.0E00;

	vector<double> n;
	n.reserve(10);
	n.push_back (0.11670521452767E04);
	n.push_back (-0.72421316703206E06);
	n.push_back (-0.17073846940092E02);
	n.push_back (0.12020824702470E05);
	n.push_back (-0.32325550322333E07);
	n.push_back (0.14915108613530E02);
	n.push_back (-0.48232657361591E04);
	n.push_back (0.40511340542057E06);
	n.push_back (-0.23855557567849E00);
	n.push_back (0.65017534844798E03);

	V = t_atm/t + (n[8]/((t_atm/t)-n[9]));
	A = pow(V,2) + n[0]*V + n[1];
	B = n[2]*pow(V,2) + n[3]*V + n[4];
	C = n[5]*pow(V,2) + n[6]*V + n[7];
	p_sat = p*pow((2*C/(-B + sqrt(pow(B,2)-4*A*C))),4);

	return p_sat;
}

vector<double> ObliczWektorX1(double &p_atm, double &p_nas, 
	double &f, vector<double> ps, vector<double> a1) {

	double p_czas, x_H2O, x_ps;
	vector<double> vect;

	f = f/100;
	p_czas = f*p_nas;
	x_H2O = p_czas/p_atm;
	x_ps = 1 - x_H2O;

	vect.push_back (ps[0]*x_ps + a1[0]*x_H2O);
	vect.push_back (ps[1]*x_ps + a1[1]*x_H2O);
	vect.push_back (ps[2]*x_ps + a1[2]*x_H2O);
	vect.push_back (ps[3]*x_ps + a1[3]*x_H2O);
	vect.push_back (ps[4]*x_ps + a1[4]*x_H2O);
	vect.push_back (ps[5]*x_ps + a1[5]*x_H2O);
	vect.push_back (ps[6]*x_ps + a1[6]*x_H2O);
	vect.push_back (ps[7]*x_ps + a1[7]*x_H2O);
	vect.push_back (ps[8]*x_ps + a1[8]*x_H2O);

	return vect;
}

Wynik funkcji ObliczCisnienieNasycenia(t_atm) jest poprawny, natomiast wektora X1 już nie. Prawidłowe wyniki wektora X1 powinny wyglądać następująco:

9.598721066535E+003
-189.025823194585E+000
4.967971055702E+000
-5.629869836106E-003
10.486675165748E-006
-7.810046822116E-009
2.149655802455E-012
-508.73338430187E+000
-3.7077387928E+000

Czyżby różnice w wynikach były spowodowane jakimś ograniczeniem STLowskiego vector'a ?

Zapomniałeś ustawić wartości p_atm, aktualnie ma jakąś losową wartość. Po ustawieniu np. na 101300.0 wyniki są prawie identycznie z tym, co podałeś.

EDIT:
[1] Dlaczego podajesz wszystkie parametry jako nie-stałe referencje (oprócz wektorów, które przekazujesz przez wartość)? Lepiej jest tak:

double ObliczCisnienieNasycenia(double t_atm);
vector<double> ObliczWektorX1(double p_atm, double p_nas, double f, const vector<double>& ps, const vector<double>& a1);

Aktualnie w funkcji ObliczWektorX1 zmieniasz przekazaną wartość f, co raczej nie powinno mieć miejsca; przekazywanie przez wartość chroni przed takimi błędami.

[2] Dlaczego w funkcji ObliczWektorX1 brak pętli? Nie lepiej zapisać tak?

vector<double> vect(9);

// ...

for (size_t i = 0; i < 9; ++i)
	vect[i] = ps[i] * x_ps + a1[i] * x_H2O;

Rzeczywiście, po ustawieniu p_atm na 101325 wyniki są identyczne. Dzięki [browar] Tylko dlaczego mimo wszystko program się wykonał ? Przecież w funkcji ObliczWektorX1() jest działanie x_H2O = p_czas/p_atm, a następnie x_H2O jest używane do obliczenia wektora vect. Przyznam że dopiero od niedawna zacząłem przygodę z C++, ale na logikę to to nie powinno przejść.

rincewind napisał(a)

[1] Dlaczego podajesz wszystkie parametry jako nie-stałe referencje (oprócz wektorów, które przekazujesz przez wartość)? Lepiej jest tak:

double ObliczCisnienieNasycenia(double t_atm);
vector<double> ObliczWektorX1(double p_atm, double p_nas, double f, const vector<double>& ps, const vector<double>& a1);

Aktualnie w funkcji ObliczWektorX1 zmieniasz przekazaną wartość f, co raczej nie powinno mieć miejsca; przekazywanie przez wartość chroni przed takimi błędami.
</cpp>

A nie lepiej przesyłać przez referencję ? Przesyłanie przez wartość powoduje przecież kopiowanie zmiennych, więc wykonuje się wolniej. Dzięki za wszelkie uwagi, naprawdę bardzo mi pomogłeś, jeśli coś byś jeszcze zmienił / rozwiązał inaczej, chętnie skorzystam z twego doświadczenia.

Nosferatu, przekazywanie przez wartość jest wolniejsze. Ale przez referencję modyfikując argument modyfikujesz oryginalną zmienną. Poza tym przez referencję nie można przesyłać stałych, chyba że w funkcji będzie zaznaczone że argument jest const i nie można go zmieniać.

Nosferatu napisał(a)

A nie lepiej przesyłać przez referencję ? Przesyłanie przez wartość powoduje przecież kopiowanie zmiennych, więc wykonuje się wolniej. Dzięki za wszelkie uwagi, naprawdę bardzo mi pomogłeś, jeśli coś byś jeszcze zmienił / rozwiązał inaczej, chętnie skorzystam z twego doświadczenia.

Owszem, lepiej przez referencję, ale nie ma potrzeby jej nadużywać. Typy podstawowe można spokojnie słać przez wartość. Zauważ, że w deklaracji zmieniłem vectory tak, aby były przesyłane przez stałą referencję, bo są nieco większymi obiektami i nie ma większego sensu ich kopiować.

Nosferatu napisał(a)

Tylko dlaczego mimo wszystko program się wykonał ? Przecież w funkcji ObliczWektorX1() jest działanie x_H2O = p_czas/p_atm, a następnie x_H2O jest używane do obliczenia wektora vect. Przyznam że dopiero od niedawna zacząłem przygodę z C++, ale na logikę to to nie powinno przejść.

To, że p_atm nie ma sensownej wartości nie znaczy, że nie ma żadnej wartości. Zmienne na stosie nie są inicjalizowane przez żadną wartość i początkowo mają wartość taką, jaka jest w pamięci w miejscu, które zostało im przydzielone. Z reguły są to śmieci, dlatego dobrą praktyką jest inicjalizacja przy deklaracji, np.:

double p_atm = 0.0,
       t_atm = 0.0,
       f = 0.0,
       p_nas = 0.0,
       x_p = 0.0;

Dzięki chłopaki za odpowiedzi. Teraz progs wygląda tak:

#include <iostream>
#include <cmath>
#include <vector>

using namespace std;

//prototypy funkcji
double ObliczCisnienieNasycenia(double t_atm);
vector<double> ObliczWektorX1(double p_atm, double p_nas, 
	double f, const vector<double>& ps, const vector<double>& a1);

int main()
{
	double p_atm, t_atm, f, p_nas, x_p;
	
	vector<double> ps;		//wektor współczynników pow. suchego
	ps.reserve(9);
	ps.push_back (1.009950160E+04);
	ps.push_back (-1.968275610E+02);
	ps.push_back (5.009155110E+00);
	ps.push_back (-5.761013730E-03);
	ps.push_back (1.066859930E-05);
	ps.push_back (-7.940297970E-09);
	ps.push_back (2.185231910E-12);
	ps.push_back (-1.767967310E+02);
	ps.push_back (-3.921504225E+00);

	vector<double> a1;		//wektor współczynników wody
	a1.reserve(9);
	a1.push_back (-3.947960830E+04);
	a1.push_back (5.755731020E+02);
	a1.push_back (9.317826530E-01);
	a1.push_back (7.222712860E-03);
	a1.push_back (-7.342557370E-06);
	a1.push_back (4.955043490E-09);
	a1.push_back (-1.336933246E-12);
	a1.push_back (-3.303974310E+04);
	a1.push_back (1.724205775E+01);

	p_atm = 101325;
	t_atm = 288.15;
	f = 60;
	f = f / 100;

	p_nas = ObliczCisnienieNasycenia(t_atm);
	cout << "Cisnienie nasycenia dla 15 *C wynosi: " << p_nas << " Pa" << endl;
	
	vector<double> X1;
	X1 = ObliczWektorX1(p_atm, p_nas, f, ps, a1);
	
	for(int i = 0; i < 9; i++) {
	cout <<  X1[i] << endl;
	}

	system("pause");
	return 0;
}

double ObliczCisnienieNasycenia(double t_atm) {

	double V, A, B, C, p_sat;
	const double p = 1.0E06, t = 1.0E00;

	vector<double> n;
	n.reserve(10);
	n.push_back (0.11670521452767E04);
	n.push_back (-0.72421316703206E06);
	n.push_back (-0.17073846940092E02);
	n.push_back (0.12020824702470E05);
	n.push_back (-0.32325550322333E07);
	n.push_back (0.14915108613530E02);
	n.push_back (-0.48232657361591E04);
	n.push_back (0.40511340542057E06);
	n.push_back (-0.23855557567849E00);
	n.push_back (0.65017534844798E03);

	V = t_atm/t + (n[8]/((t_atm/t)-n[9]));
	A = pow(V,2) + n[0]*V + n[1];
	B = n[2]*pow(V,2) + n[3]*V + n[4];
	C = n[5]*pow(V,2) + n[6]*V + n[7];
	p_sat = p*pow((2*C/(-B + sqrt(pow(B,2)-4*A*C))),4);

	return p_sat;
}

vector<double> ObliczWektorX1(double p_atm, double p_nas, 
	double f, const vector<double>& ps, const vector<double>& a1) {

	double p_czas, x_H2O, x_ps;
	vector<double> vect(9);

	p_czas = f*p_nas;
	x_H2O = p_czas/p_atm;
	x_ps = 1 - x_H2O;

	for(size_t i = 0; i < 9; ++i)
		vect[i] = ps[i]*x_ps + a1[i]*x_H2O;

	return vect;
}

Wstawiłem f = f/100 do main'a. Mam nadzieję że tak się przesyła przez wartość :) Aha, nie wiem po co jest to size_t. Czy teraz jest wszystko w porządku ?

Teraz wszystko OK. size_t jest "typem rozmiaru" -- właściwie to jest liczba całkowita bez znaku. Nie ma przeszkód, żeby zamiast niego wstawić int. W Visual C++ 2008 size_t jest zdefiniowane następująco:

typedef unsigned int size_t;

Użyłem size_t z przyzwyczajenia, choć można użyć jeszcze bardziej "pure STL" metody i zamiast size_t dać vector<double>::size_type. :P Ale int jest w zupełności wystarczający do większości zastosowań.

Nie, nie może być int. Nigdy nie stosuj rozmiaru/indeksu mogącego przyjmować wartości ujemne.

@deus, to sprawa dość dyskusyjna. Właściwie powinno się stosować taki sam typ, jakim indeksowany jest wektor. Akurat std::vector ma indeksy unsigned, więc takiego samego typu licznik powinno się stosować -- mój błąd.

Ale indeksy unsigned też mogą trochę namieszać i trzeba uważać. Na przykład taka pętla:

for (unsigned int i = 10; i >= 0; --i)
   // do sth

I mamy pętlę nieskończoną. Powiedziałbym wręcz, że indeksy signed są bezpieczniejszym wyjściem, ale że STL korzysta z unsigned, nie powinno się ich mieszać z signed, bo to prowadzi do jeszcze gorszych błędów.

Mówiłem o indeksach tablic/wtf. Użycie typów signed chroni przed idiotycznymi błędami, które i tak natychmiast wyłapiesz z użyciem debuggera w wypadku wartości unsigned. W innych wypadkach jest to zwyczajnie groźne, niebezpieczne bo może tworzyć w programie podatności.

Czyli najlepiej będzie zostawić tak jak jest (size_t) czy zamiast tego dać vector<double>::size_type ? Rozumiem że int'a tam nie powinienem wstawiać.

To jest w zasadzie obojętne, ważne, żeby był to typ całkowity bez znaku. Może być size_t, vector<double>::size_type, unsigned int... Do wyboru do koloru. :)