|
|
★ ★ ★ ★ ★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
resonant(金币+5):欢迎交流，感谢回帖:-) 2010-04-21 18:34
liangzidou(金币+5):谢谢你的热情帮助，因为你回帖时金币已经赠送完了，就给你追加5个，币不多，略表谢意。 2010-04-22 18:00
wangen994(程序强帖+1):感谢你的参与，呵呵，欢迎常来 2010-05-09 21:49:52
看到你说热力学，来个C/C++版本的，如果你能知道压力对压缩因子/体积/密度的导数（二阶导数）信息的话，用迭代法求也不错，毕竟有时候三角函数计算也很耗时
#include
#include
#include
using namespace std;
typedef unsigned char Boolean;
#define False ((Boolean)'\0')
#define True ((Boolean)'\1')
typedef int Index;
typedef double Real;
// The last parameter denotes complex roots is allowed or not
// X^3 + Coeff[0]*X^2 + Coeff[1]*X + Coeff[2] = 0
Index CubicRoot(Real *Coeff, Real *Root, Boolean Complex=True)
{
Real p=Coeff[0], q=Coeff[1], r=Coeff[2];
Real Den_3, A, A3, B, B2, D, sqrt_D, M, N;
Den_3 = 0.33333333333333333333333333;
Real Den_27 = Den_3*Den_3*Den_3;
A = Den_3 * (3.0*q - p*p);
A3 = Den_27*A*A*A;
B = 2.0*Den_27*p*p*p - Den_3*p*q + r;
B2 = 0.25*B*B;
D = A3 + B2;
if (D >= 0.0)
{
sqrt_D = sqrt(D);
M = pow((-0.5*B + sqrt_D), Den_3);
N = pow((-0.5*B - sqrt_D), Den_3);
if (D > 0.0) // 1 root
{
Root[0] = M + N;
nRoot = 1;
if (bComplex)
{
Root[1] = -0.5*(M + N);
Root[2] = 0.5*sqrt(3.0)*(M - N);
nRoot = 3;
}
}
else // 2 roots
{
Root[0] = M + N;
Root[1] = Root[2] = -0.5*Root[0];
nRoot = 2;
}
}
else // 3 real distinct roots
{
Real deg = sqrt(-B2/A3);
Real PI_2 = 1.5707963267949; // 90C
Real PI_3_2 = 2.0943951023932; // 120C
Real Theta;
if (B >= 0.0)
{
Theta = (PI_2 + atan(deg/sqrt(1.0 - deg*deg)))*Den_3;
}
else
{
Theta = atan(sqrt(1.0 - deg*deg)/deg)*Den_3;
}
Real factor = 2.0*sqrt(-Den_3*A);
for (Index i=0; i<3; i++)
{
Root = factor*cos(Theta + PI_3_2*i);
}
nRoot = 3;
}
for (Index i=0; i<3; i++)
{
Root -= Den_3*p;
}
return nRoot;
}
int main(int argc, char* argv[])
{
string Holding;
Real a[3] = {-0.9964370251388, 0.0345649889308, -0.000135939438932387};
Real root[3];
Index nRoot = CubicRoot(a, root);
cout<<"Root 1: "<
cout<<"Root 2: "<
cout<<"Root 3: "<
cin>>Holding;
return 0;
} |
|
回复此楼