|
|
|
@ -762,6 +762,228 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests |
|
|
|
return Broyden.FindRoot(fw, initialGuess, accuracy, maxIterations)[0]; |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Twoeq1() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Twoeq1.htm
|
|
|
|
// not solvable with this method
|
|
|
|
Func<double[], double[]> fa1 = x => |
|
|
|
{ |
|
|
|
double D = x[0]; |
|
|
|
double fF = x[1]; |
|
|
|
const double dp = 103000; |
|
|
|
const double L = 100; |
|
|
|
const double T = 25 + 273.15; |
|
|
|
const double Q = 0.0025; |
|
|
|
const double pi = 3.1416; |
|
|
|
const double rho = 46.048 + T * (9.418 + T * (-0.0329 + T * (4.882e-5 - T * 2.895e-8))); |
|
|
|
double vis = Math.Exp(-10.547 + 541.69 / (T - 144.53)); |
|
|
|
double v = Q / (pi * D * D / 4); |
|
|
|
double kvis = vis / rho; |
|
|
|
double Re = v * D / kvis; |
|
|
|
double fD = -dp / rho + 2 * fF * v * v * L / D; |
|
|
|
double ffF = (Re < 2100) ? (fF - 16 / Re) : (fF - 1 / Math.Pow(4 * Math.Log(Re * Math.Sqrt(fF)) - 0.4, 2)); |
|
|
|
return new double[2] { fD, ffF }; |
|
|
|
}; |
|
|
|
|
|
|
|
Assert.Throws<NonConvergenceException>(() => Broyden.FindRoot(fa1, new double[2] { 0.04, 0.001 }, 1e-14)); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Twoeq2() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Twoeq2.htm
|
|
|
|
// with this method this test is only solvable with reduced accuracy
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double x = xa[0]; |
|
|
|
double T = xa[1]; |
|
|
|
double k = 0.12 * Math.Exp(12581 * (T - 298) / (298 * T)); |
|
|
|
double fx = 120 * x - 75 * k * (1 - x); |
|
|
|
double fT = -x * (873 - T) + 11.0 * (T - 300); |
|
|
|
return new double[2] { fx, fT }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[2] { 1, 400 }, 1e-6); |
|
|
|
Assert.AreEqual(0.9638680512795, r[0], 1e-5); |
|
|
|
Assert.AreEqual(346.16369814640, r[1], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-6); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-7); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Twoeq3() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Twoeq3.htm
|
|
|
|
// with this method this test is only solvable with slightly reduced accuracy
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double x = xa[0]; |
|
|
|
double T = xa[1]; |
|
|
|
double k = Math.Exp(-149750 / T + 92.5); |
|
|
|
double Kp = Math.Exp(42300 / T - 24.2 + 0.17 * Math.Log(T)); |
|
|
|
double fx = k * Math.Sqrt(1 - x) * ((0.91 - 0.5 * x) / (9.1 - 0.5 * x) - x * x / ((1 - x) * (1 - x) * Kp)); |
|
|
|
double fT = T * (1.84 * x + 77.3) - 43260 * x - 105128; |
|
|
|
return new double[2] { fx, fT }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[2] { 0.5, 1700 }, 1e-11); |
|
|
|
Assert.AreEqual(0.5333728995523, r[0], 1e-5); |
|
|
|
Assert.AreEqual(1637.70322946500, r[1], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-11); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-14); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Twoeq4a() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Twoeq4a.htm
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double G21 = xa[0]; |
|
|
|
double G12 = xa[1]; |
|
|
|
const double t = 58.7; |
|
|
|
const double pw1 = 21; |
|
|
|
const double x1 = (pw1 / 46.07) / (pw1 / 46.07 + (100 - pw1) / 86.18); |
|
|
|
double P2 = Math.Pow(10, 6.87776 - 1171.53 / (224.366 + t)); |
|
|
|
double P1 = Math.Pow(10, 8.04494 - 1554.3 / (222.65 + t)); |
|
|
|
double gamma2 = 760 / P2; |
|
|
|
double gamma1 = 760 / P1; |
|
|
|
const double x2 = 1 - x1; |
|
|
|
double t1 = x1 + x2 * G12; |
|
|
|
double t2 = x2 + x1 * G21; |
|
|
|
double g2calc = Math.Exp(-Math.Log(t2) - x1 * (G12 * t2 - G21 * t1) / (t1 * t2)); |
|
|
|
double g1calc = Math.Exp(-Math.Log(t1) + x2 * (G12 * t2 - G21 * t1) / (t1 * t2)); |
|
|
|
double fG21 = Math.Log(gamma2) + Math.Log(t2) + x1 * (G12 * t2 - G21 * t1) / (t1 * t2); |
|
|
|
double fG12 = Math.Log(gamma1) + Math.Log(t1) - x2 * (G12 * t2 - G21 * t1) / (t1 * t2); |
|
|
|
return new double[2] { fG21, fG12 }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[2] { 0.5, 0.5 }, 1e-14); |
|
|
|
Assert.AreEqual(0.3017535930592, r[0], 1e-5); |
|
|
|
Assert.AreEqual(0.0785888476379, r[1], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-14); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Twoeq4b() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Twoeq4b.htm
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double G21 = xa[0]; |
|
|
|
double G12 = xa[1]; |
|
|
|
const double t = 58.7; |
|
|
|
const double pw1 = 21; |
|
|
|
const double x1 = (pw1 / 46.07) / (pw1 / 46.07 + (100 - pw1) / 86.18); |
|
|
|
double P2 = Math.Pow(10, 6.87776 - 1171.53 / (224.366 + t)); |
|
|
|
double P1 = Math.Pow(10, 8.04494 - 1554.3 / (222.65 + t)); |
|
|
|
double gamma2 = 760 / P2; |
|
|
|
double gamma1 = 760 / P1; |
|
|
|
const double x2 = 1 - x1; |
|
|
|
double t1 = x1 + x2 * G12; |
|
|
|
double t2 = x2 + x1 * G21; |
|
|
|
double g2calc = Math.Exp(-Math.Log(t2) - x1 * (G12 * t2 - G21 * t1) / (t1 * t2)); |
|
|
|
double g1calc = Math.Exp(-Math.Log(t1) + x2 * (G12 * t2 - G21 * t1) / (t1 * t2)); |
|
|
|
double fG21 = t1 * t2 * (Math.Log(gamma2) + Math.Log(t2)) + x1 * (G12 * t2 - G21 * t1); |
|
|
|
double fG12 = t1 * t2 * (Math.Log(gamma1) + Math.Log(t1)) - x2 * (G12 * t2 - G21 * t1); |
|
|
|
return new double[2] { fG21, fG12 }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[2] { 0.8, 0.8 }, 1e-14); |
|
|
|
Assert.AreEqual(0.3017535528355, r[0], 1e-5); |
|
|
|
Assert.AreEqual(0.0785888934885, r[1], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-14); |
|
|
|
} |
|
|
|
|
|
|
|
/* |
|
|
|
* Test converges, but to a point at A=1.7455, B=2.59. Error in test case? |
|
|
|
[Test] |
|
|
|
public void Twoeq5a() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Twoeq5a.htm
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double A = xa[0]; |
|
|
|
double B = xa[1]; |
|
|
|
const double t = 70.9; |
|
|
|
const double pw1 = 49; |
|
|
|
const double x1 = (pw1 / 46.07) / (pw1 / 46.07 + (100 - pw1) / 100.2); |
|
|
|
double P2 = Math.Pow(10, 6.9024 - 1268.115 / (216.9 + t)); |
|
|
|
double P1 = Math.Pow(10, 8.04494 - 1554.3 / (222.65 + t)); |
|
|
|
double gamma1 = 760 / P1; |
|
|
|
const double x2 = 1 - x1; |
|
|
|
double gamma2 = 760 / P2; |
|
|
|
double g1calc = Math.Pow(10, A * Math.Pow(x2 / (A * x1 / B + x2), 2)); |
|
|
|
double g2calc = Math.Pow(10, B * Math.Pow(x1 / (x1 + B * x2 / A), 2)); |
|
|
|
double fA = Math.Log(gamma1) - A * Math.Pow(x2 / (A * x1 / B + x2), 2); |
|
|
|
double fB = Math.Log(gamma2) - B * Math.Pow(x1 / (x1 + B * x2 / A), 2); |
|
|
|
return new double[2] { fA, fB }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[2] { 1.0, 1.0 }, 1e-14); |
|
|
|
Assert.AreEqual(0.7580768059470, r[0], 1e-5); |
|
|
|
Assert.AreEqual(1.1249034445330, r[1], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-14); |
|
|
|
} |
|
|
|
*/ |
|
|
|
|
|
|
|
/* |
|
|
|
* Test converges, but to a point at A=1.7455, B=2.59. Error in test case? |
|
|
|
[Test] |
|
|
|
public void Twoeq5b() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Twoeq5b.htm
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double A = xa[0]; |
|
|
|
double B = xa[1]; |
|
|
|
const double t = 70.9; |
|
|
|
const double pw1 = 49; |
|
|
|
const double x1 = (pw1 / 46.07) / (pw1 / 46.07 + (100 - pw1) / 100.2); |
|
|
|
double P2 = Math.Pow(10, 6.9024 - 1268.115 / (216.9 + t)); |
|
|
|
double P1 = Math.Pow(10, 8.04494 - 1554.3 / (222.65 + t)); |
|
|
|
double gamma1 = 760 / P1; |
|
|
|
const double x2 = 1 - x1; |
|
|
|
double gamma2 = 760 / P2; |
|
|
|
double g1calc = Math.Pow(10, A * Math.Pow(x2 / (A * x1 / B + x2), 2)); |
|
|
|
double g2calc = Math.Pow(10, B * Math.Pow(x1 / (x1 + B * x2 / A), 2)); |
|
|
|
double fA = Math.Log(gamma1) * Math.Pow(A * x1 / B + x2, 2) - A * Math.Pow(x2, 2); |
|
|
|
double fB = Math.Log(gamma2) * Math.Pow(x1 + B * x2 / A, 2) - B * Math.Pow(x1, 2); |
|
|
|
return new double[2] { fA, fB }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[2] { 1.0, 1.0 }, 1e-14); |
|
|
|
Assert.AreEqual(0.7580768919420, r[0], 1e-5); |
|
|
|
Assert.AreEqual(1.1249035089540, r[1], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-14); |
|
|
|
} |
|
|
|
*/ |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Twoeq6() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Twoeq7.htm
|
|
|
|
Func<double[], double[]> fa1 = x => |
|
|
|
{ |
|
|
|
double x1 = x[0]; |
|
|
|
double x2 = x[1]; |
|
|
|
double f1 = x1 / (1 - x1) - 5 * Math.Log(0.4 * (1 - x1) / x2) + 4.45977; |
|
|
|
double f2 = x2 - (0.4 - 0.5 * x1); |
|
|
|
return new double[2] { f1, f2 }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[2] { 0.9, 0.5 }, 1e-14); |
|
|
|
Assert.AreEqual(0.7573962468236, r[0], 1e-5); |
|
|
|
Assert.AreEqual(0.0213018765882, r[1], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-14); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Twoeq7() |
|
|
|
{ |
|
|
|
@ -786,5 +1008,370 @@ namespace MathNet.Numerics.UnitTests.RootFindingTests |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-14); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Twoeq8() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Twoeq8.htm
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double rp = xa[0]; |
|
|
|
double x = xa[1]; |
|
|
|
double frp = rp - 0.327 * Math.Pow(x, 0.804) * Math.Exp(-5230 / (1.987 * (373 + 1.84e6 * rp))); |
|
|
|
double fx = x - (0.06 - 161 * rp); |
|
|
|
return new double[2] { frp, fx }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[2] { 0.0001, 0.01 }, 1e-14); |
|
|
|
Assert.AreEqual(0.0003406054400, r[0], 1e-5); |
|
|
|
Assert.AreEqual(0.0051625241669, r[1], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-14); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Twoeq9() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Twoeq9.htm
|
|
|
|
// not solvable with this method
|
|
|
|
Func<double[], double[]> fa1 = x => |
|
|
|
{ |
|
|
|
double fF = x[0]; |
|
|
|
double u = x[1]; |
|
|
|
const double L = 6000; |
|
|
|
const double D = 0.505; |
|
|
|
const double rho = 53; |
|
|
|
const double g = 32.2; |
|
|
|
const double eps = .00015; |
|
|
|
double Re = rho * D * u / (13.2 * 0.000672); |
|
|
|
double ffF = fF - 1 / Math.Pow(2.28 - 4 * Math.Log(eps / D + 4.67 / (Re * Math.Sqrt(fF))), 2); |
|
|
|
double fu = 133.7 - (2 * fF * rho * u * u * L / D + rho * g * 200) / (g * 144); |
|
|
|
return new double[2] { ffF, fu }; |
|
|
|
}; |
|
|
|
|
|
|
|
Assert.Throws<NonConvergenceException>(() => Broyden.FindRoot(fa1, new double[2] { 0.1, 10.0 }, 1e-14)); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Twoeq10() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Twoeq10.htm
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double t12 = xa[0]; |
|
|
|
double t21 = xa[1]; |
|
|
|
const double alpha = 0.4; |
|
|
|
double c1 = Math.Exp(-alpha * t21); |
|
|
|
double c2 = Math.Exp(-2 * alpha * t21); |
|
|
|
double c3 = Math.Exp(-alpha * t12); |
|
|
|
double c4 = Math.Exp(-2 * alpha * t12); |
|
|
|
const double x1 = 0.5; |
|
|
|
const double x2 = 1 - x1; |
|
|
|
|
|
|
|
double ft12 = 1 / (x1 * x2) - 2 * t21 * c2 / Math.Pow(x1 + x2 * c1, 3) - 2 * t12 * c4 / Math.Pow(x2 + x1 * c3, 3); |
|
|
|
double ft21 = (x1 - x2) / Math.Pow(x1 * x2, 2) + 6 * t21 * c2 * (1 - c1) / Math.Pow(x1 + x2 * c1, 4) + 6 * t12 * c4 * (c3 - 1) / Math.Pow(x2 + x1 * c3, 4); |
|
|
|
return new double[2] { ft12, ft21 }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[2] { 0.1, 0.1 }, 1e-14); |
|
|
|
Assert.AreEqual(1.6043843214350, r[0], 1e-5); |
|
|
|
Assert.AreEqual(1.6043843214350, r[1], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-14); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Threeq1() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Threeq1.htm
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double t = xa[0]; |
|
|
|
double x1 = xa[1]; |
|
|
|
double x2 = xa[2]; |
|
|
|
const double y2 = 0.8; |
|
|
|
const double y1 = 0.2; |
|
|
|
double p1 = Math.Pow(10, 7.62231 - 1417.9 / (191.15 + t)); |
|
|
|
double p2 = Math.Pow(10, 8.10765 - 1750.29 / (235 + t)); |
|
|
|
const double B = 0.7; |
|
|
|
const double A = 1.7; |
|
|
|
double gamma2 = Math.Pow(10, B * x1 * x1 / Math.Pow(x1 + B * x2 / A, 2)); |
|
|
|
double gamma1 = Math.Pow(10, A * x2 * x2 / Math.Pow(A * x1 / B + x2, 2)); |
|
|
|
double k2 = gamma2 * p2 / 760; |
|
|
|
double k1 = gamma1 * p1 / 760; |
|
|
|
double ft = x1 + x2 - 1; |
|
|
|
double fx1 = x1 - y1 / k1; |
|
|
|
double fx2 = x2 - y2 / k2; |
|
|
|
return new double[3] { ft, fx1, fx2 }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[3] { 100, 0.2, 0.8 }, 1e-14); |
|
|
|
Assert.AreEqual(93.96706523770, r[0], 1e-5); |
|
|
|
Assert.AreEqual(0.0078754574659, r[1], 1e-5); |
|
|
|
Assert.AreEqual(0.9921245425339, r[2], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[2], 1e-14); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Threeq2() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Threeq2.htm
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double x1 = xa[0]; |
|
|
|
double x2 = xa[1]; |
|
|
|
double alpha = xa[2]; |
|
|
|
const double t = 88.538; |
|
|
|
const double B = 0.7; |
|
|
|
const double A = 1.7; |
|
|
|
const double z1 = 0.2; |
|
|
|
const double z2 = 0.8; |
|
|
|
double p1 = Math.Pow(10, 7.62231 - 1417.9 / (191.15 + t)); |
|
|
|
double p2 = Math.Pow(10, 8.10765 - 1750.29 / (235 + t)); |
|
|
|
double gamma2 = Math.Pow(10, B * x1 * x1 / Math.Pow(x1 + B * x2 / A, 2)); |
|
|
|
double gamma1 = Math.Pow(10, A * x2 * x2 / Math.Pow(A * x1 / B + x2, 2)); |
|
|
|
double k1 = gamma1 * p1 / 760; |
|
|
|
double k2 = gamma2 * p2 / 760; |
|
|
|
double y1 = k1 * x1; |
|
|
|
double y2 = k2 * x2; |
|
|
|
double fx1 = x1 - z1 / (1 + alpha * (k1 - 1)); |
|
|
|
double fx2 = x2 - z2 / (1 + alpha * (k2 - 1)); |
|
|
|
double falpha = x1 + x2 - (y1 + y2); |
|
|
|
return new double[3] { fx1, fx2, falpha }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[3] { 0, 1, 0.5 }, 1e-14); |
|
|
|
Assert.AreEqual(0.0226974766367, r[0], 1e-5); |
|
|
|
Assert.AreEqual(0.9773025233633, r[1], 1e-5); |
|
|
|
Assert.AreEqual(0.5322677863643, r[2], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[2], 1e-14); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Threeq3() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Threeq3.htm
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double T = xa[0]; |
|
|
|
double Ca = xa[1]; |
|
|
|
double Tj = xa[2]; |
|
|
|
double k = 7.08e10 * Math.Exp(-30000 / 1.9872 / T); |
|
|
|
const double rhocp = 50 * 0.75; |
|
|
|
const double F = 40; |
|
|
|
const double T0 = 530; |
|
|
|
const double V = 48; |
|
|
|
const double U = 150; |
|
|
|
const double A = 250; |
|
|
|
const double Ca0 = 0.55; |
|
|
|
const double Fj = 49.9; |
|
|
|
const double Tj0 = 530; |
|
|
|
double fT = F * (T0 - T) / V + 30000 * k * Ca / rhocp - U * A * (T - Tj) / (rhocp * V); |
|
|
|
double fCa = F * (Ca0 - Ca) / V - k * Ca; |
|
|
|
double fTj = Fj * (Tj0 - Tj) / 3.85 + U * A * (T - Tj) / (62.3 * 1.0 * 3.85); |
|
|
|
return new double[3] { fT, fCa, fTj }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[3] { 600, 0.1, 600 }, 1e-11); |
|
|
|
Assert.AreEqual(590.34979512380, r[0], 1e-5); |
|
|
|
Assert.AreEqual(0.3301868979161, r[1], 1e-5); |
|
|
|
Assert.AreEqual(585.72976766210, r[2], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-12); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[2], 1e-11); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Threeq4a() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Threeq4a.htm
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double CD = xa[0]; |
|
|
|
double CX = xa[1]; |
|
|
|
double CZ = xa[2]; |
|
|
|
double CY = CX + CZ; |
|
|
|
double CC = CD - CY; |
|
|
|
const double CA0 = 1.5; |
|
|
|
const double CB0 = 1.5; |
|
|
|
double CA = CA0 - CD - CZ; |
|
|
|
double CB = CB0 - CD - CY; |
|
|
|
const double KC1 = 1.06; |
|
|
|
const double KC2 = 2.63; |
|
|
|
const double KC3 = 5; |
|
|
|
double fCD = CC * CD / (CA * CB) - KC1; |
|
|
|
double fCX = CX * CY / (CB * CC) - KC2; |
|
|
|
double fCZ = CZ / (CA * CX) - KC3; |
|
|
|
return new double[3] { fCD, fCX, fCZ }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[3] { 0.7, 0.2, 0.4 }, 1e-14); |
|
|
|
Assert.AreEqual(0.7053344059695, r[0], 1e-5); |
|
|
|
Assert.AreEqual(0.1777924200537, r[1], 1e-5); |
|
|
|
Assert.AreEqual(0.3739765850146, r[2], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[2], 1e-14); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Threeq4b() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Threeq4b.htm
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double CD = xa[0]; |
|
|
|
double CX = xa[1]; |
|
|
|
double CZ = xa[2]; |
|
|
|
double CY = CX + CZ; |
|
|
|
double CC = CD - CY; |
|
|
|
const double CA0 = 1.5; |
|
|
|
const double CB0 = 1.5; |
|
|
|
double CA = CA0 - CD - CZ; |
|
|
|
double CB = CB0 - CD - CY; |
|
|
|
const double KC1 = 1.06; |
|
|
|
const double KC2 = 2.63; |
|
|
|
const double KC3 = 5; |
|
|
|
double fCD = CC * CD - KC1 * (CA * CB); |
|
|
|
double fCX = CX * CY - KC2 * (CB * CC); |
|
|
|
double fCZ = CZ - KC3 * (CA * CX); |
|
|
|
return new double[3] { fCD, fCX, fCZ }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[3] { 0.7, 0.2, 0.4 }, 1e-14); |
|
|
|
Assert.AreEqual(0.7053344059695, r[0], 1e-5); |
|
|
|
Assert.AreEqual(0.1777924200537, r[1], 1e-5); |
|
|
|
Assert.AreEqual(0.3739765850146, r[2], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[2], 1e-14); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Threeq5() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Threeq5.htm
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double X = xa[0]; |
|
|
|
double T = xa[1]; |
|
|
|
double h = xa[2]; |
|
|
|
double k1 = 3e5 * Math.Exp(-5000 / T); |
|
|
|
double k2 = 6e7 * Math.Exp(-7500 / T); |
|
|
|
const double F0 = 1; |
|
|
|
const double T0 = 300; |
|
|
|
double fX = -0.16 * X * F0 / h + k1 * (1 - X) - k2 * X; |
|
|
|
double fT = 0.16 * F0 * T0 / h - 0.16 * T * F0 / h + 5 * (k1 * (1 - X) - k2 * X); |
|
|
|
double fh = 0.16 * F0 - 0.4 * Math.Sqrt(h); |
|
|
|
return new double[3] { fX, fT, fh }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[3] { 0.5, 500, 0.5 }, 1e-13); |
|
|
|
Assert.AreEqual(0.0171035426725, r[0], 1e-5); |
|
|
|
Assert.AreEqual(300.08551771340, r[1], 1e-5); |
|
|
|
Assert.AreEqual(0.16, r[2], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-13); |
|
|
|
Assert.AreEqual(0, fa1(r)[2], 1e-14); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Threeq6() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Threeq6.htm
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double CA = xa[0]; |
|
|
|
double CB = xa[1]; |
|
|
|
double T = xa[2]; |
|
|
|
double k1 = 11 * Math.Exp(-4180 / (8.314 * (T + 273.16))); |
|
|
|
double k2 = 172.2 * Math.Exp(-34833 / (8.314 * (T + 273.16))); |
|
|
|
const double Q = 5.1E6; |
|
|
|
double fCA = 0.1 * (1 - CA) - k1 * CA * CA; |
|
|
|
double fCB = -0.1 * CB + k1 * CA * CA - k2 * CB; |
|
|
|
double fT = 0.1 * (25 - T) - 418 * k1 * CA * CA - 418 * k2 * CB + Q * 1e-5; |
|
|
|
return new double[3] { fCA, fCB, fT }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[3] { 0.5, 0.5, 500 }, 1e-13); |
|
|
|
Assert.AreEqual(0.1578109142617, r[0], 1e-5); |
|
|
|
Assert.AreEqual(0.77071354919, r[1], 1e-5); |
|
|
|
Assert.AreEqual(153.09, r[2], 1e-2); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[2], 1e-13); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Threeq7() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Threeq7.htm
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double CA = xa[0]; |
|
|
|
double CB = xa[1]; |
|
|
|
double T = xa[2]; |
|
|
|
const double T0 = 298; |
|
|
|
const double CA0 = 3; |
|
|
|
const double CB0 = 0; |
|
|
|
const double CC0 = 0; |
|
|
|
const double theta = 300; |
|
|
|
double k1 = 4e6 * Math.Exp(-60000 / (8.314 * T)); |
|
|
|
double KA = 17 * Math.Exp(-7000 / (8.314 * T)); |
|
|
|
double k2 = 3e4 * Math.Exp(-80000 / (8.314 * T)); |
|
|
|
double k2p = 3e4 * Math.Exp(-90000 / (8.314 * T)); |
|
|
|
double CC = (CC0 + theta * k2 * CB) / (1 + theta * k2p); |
|
|
|
double fCA = CA0 - CA - theta * k1 * CA / (1 + KA * CB); |
|
|
|
double fCB = CB - CB0 - (theta * k1 * CA / (1 + KA * CB) - theta * k2 * CB + theta * k2p * CC); |
|
|
|
double fT = 85 * (T - T0) + 0.02 * (T * T - T0 * T0) - ((16000 + 3 * T - 0.002 * T * T) * ((CA0 - CA) / CA0) + (30000 + 4 * T - 0.003 * T * T) * CC / CA0); |
|
|
|
return new double[3] { fCA, fCB, fT }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[3] { 3, 0, 300 }, 1e-13); |
|
|
|
Assert.AreEqual(2.7873203092940, r[0], 1e-5); |
|
|
|
Assert.AreEqual(0.2126796260201, r[1], 1e-5); |
|
|
|
Assert.AreEqual(310.212556340, r[2], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-14); |
|
|
|
Assert.AreEqual(0, fa1(r)[2], 1e-14); |
|
|
|
} |
|
|
|
|
|
|
|
[Test] |
|
|
|
public void Threeq8() |
|
|
|
{ |
|
|
|
// Test case from http://www.polymath-software.com/library/nle/Threeq8.htm
|
|
|
|
Func<double[], double[]> fa1 = xa => |
|
|
|
{ |
|
|
|
double p4 = xa[0]; |
|
|
|
double Q24 = xa[1]; |
|
|
|
double Q34 = xa[2]; |
|
|
|
double p2 = 156.6 - 0.00752 * Q24 * Q24; |
|
|
|
double p3 = 117.1 - 0.00427 * Q34 * Q34; |
|
|
|
const double D34 = 2.067 / 12; |
|
|
|
const double D24 = 1.278 / 12; |
|
|
|
const double D45 = 2.469 / 12; |
|
|
|
const double fF = 0.015; |
|
|
|
const double pi = 3.1416; |
|
|
|
double k24 = 2 * fF * 125 / (Math.Pow((60 * 7.48) * (pi * D24 * D24 / 4), 2) * D24); |
|
|
|
double k34 = 2 * fF * 125 / (Math.Pow((60 * 7.48) * (pi * D34 * D34 / 4), 2) * D34); |
|
|
|
double k45 = 2 * fF * 145 / (Math.Pow((60 * 7.48) * (pi * D45 * D45 / 4), 2) * D45); |
|
|
|
double fp4 = 70 * 32.3 - p4 * (144 * 32.2 / 62.35) + k45 * Math.Pow(Q24 + Q34, 2); |
|
|
|
double fQ24 = (p4 - p2) * (144 * 32.2 / 62.35) + k24 * Q24 * Q24; |
|
|
|
double fQ34 = (p4 - p3) * (144 * 32.2 / 62.35) + k34 * Q34 * Q34; |
|
|
|
return new double[3] { fp4, fQ24, fQ34 }; |
|
|
|
}; |
|
|
|
|
|
|
|
double[] r = Broyden.FindRoot(fa1, new double[3] { 50, 100, 100 }, 1e-11); |
|
|
|
Assert.AreEqual(57.12556038475, r[0], 1e-5); |
|
|
|
Assert.AreEqual(51.75154563498, r[1], 1e-5); |
|
|
|
Assert.AreEqual(92.91811138918, r[2], 1e-5); |
|
|
|
Assert.AreEqual(0, fa1(r)[0], 1e-12); |
|
|
|
Assert.AreEqual(0, fa1(r)[1], 1e-12); |
|
|
|
Assert.AreEqual(0, fa1(r)[2], 1e-12); |
|
|
|
} |
|
|
|
} |
|
|
|
} |
|
|
|
|
taschna
13 years ago