diyifan001 发表于 2010-4-16 16:46:47

以前无聊的时候搞得个东东（二）

接下来就是要考虑求压缩因子的问题了因为立方型状态方程为三次的，所以要考虑解三次方程的问题，传统解法为迭代法，矩阵分解法等，太复杂了，网上有个文章介绍了个公式，叫盛金公式，具体原理我搞不太懂，但通过在下的反复试验，得出的解和matlab的一样准确    作者好像叫范盛金，这年头搞基础学科的太少了，个人感觉他是个牛人。
PR方程整理成压缩因子的形式为三次方程   解方程   就可以得出系统的压缩因子

利用盛金公式求方程实根的代码：
public class commenGJ
    {
      //利用盛金公式，计算三次方程的最小实数根，感谢盛金同学
      public static double minSHJ(double a,double b,double c,double d)
      {
            double A = b*b-3*a*c;
            double B = b * c - 9 * a * d;
            double C = c * c - 3 * b * d;
            double dt = B * B - 4 * A * C;
            if (A == B && A == 0d)
            {
                return -c / b;
            }
            else if (dt > 0)
            {
                double Y1 = A * b + 3 * a * (-B + Math.Pow(dt, 0.5d)) / 2;
                double Y2 = A * b + 3 * a * (-B - Math.Pow(dt, 0.5d)) / 2;
                return (-b - Math.Sign(Y1) * Math.Pow(Math.Abs(Y1), 1d / 3) - Math.Sign(Y2) * Math.Pow(Math.Abs(Y2), 1d / 3)) / (3 * a);
            }
            else if (dt == 0)
            {
                double K = B / A;
                return Math.Min(-b / a + K, -K / 2);
            }
            else if (dt < 0)
            {
                double T = (2 * A * b - 3 * a * B) / (2 * Math.Sqrt(A * A * A));
                double st = Math.Acos(T);
                double temp1 = (-b - 2 * Math.Sqrt(A) * Math.Cos(st / 3)) / (3 * a);
                double temp2 = Math.Min((-b + Math.Sqrt(A) * (Math.Cos(st / 3) + Math.Sqrt(3) * Math.Sin(st / 3))) / (3 * a), (-b + Math.Sqrt(A) * (Math.Cos(st / 3) - Math.Sqrt(3) * Math.Sin(st / 3))) / (3 * a));
                return Math.Min(temp1, temp2);
               
            }
            else
                return 0d;
      }
      //利用盛金公式，计算三次方程的最大实数根，感谢盛金同学
      public static double maxSHJ(double a, double b, double c, double d)
      {
            double A = b * b - 3 * a * c;
            double B = b * c - 9 * a * d;
            double C = c * c - 3 * b * d;
            double dt = B * B - 4 * A * C;
            if (A == B && A == 0d)
            {
                return -c / b;
            }
            else if (dt > 0)
            {
                double Y1 = A * b + 3 * a * (-B + Math.Pow(dt, 0.5d)) / 2;
                double Y2 = A * b + 3 * a * (-B - Math.Pow(dt, 0.5d)) / 2;
                return (-b - Math.Sign(Y1) * Math.Pow(Math.Abs(Y1), 1d / 3) - Math.Sign(Y2) * Math.Pow(Math.Abs(Y2), 1d / 3)) / (3 * a);
            }
            else if (dt == 0)
            {
                double K = B / A;
                return Math.Max(-b / a + K, -K / 2);
            }
            else if (dt < 0)
            {
                double T = (2 * A * b - 3 * a * B) / (2 * Math.Sqrt(A * A * A));
                double st = Math.Acos(T);
                double temp1 = (-b - 2 * Math.Sqrt(A) * Math.Cos(st / 3)) / (3 * a);
                double temp2 = Math.Max((-b + Math.Sqrt(A) * (Math.Cos(st / 3) + Math.Sqrt(3) * Math.Sin(st / 3))) / (3 * a), (-b + Math.Sqrt(A) * (Math.Cos(st / 3) - Math.Sqrt(3) * Math.Sin(st / 3))) / (3 * a));
                return Math.Max(temp1, temp2);
            }
            else
                return 0d;
      }

      
    }

写的太仓促，欢迎大家一起讨论 参与
说不准将来远程合作开发个炮网专用的计算软件。
哈哈

zhuniuBB 发表于 2010-4-16 16:46:47

好有趣的编程。

wmm598 发表于 2010-4-16 16:46:47


这东西我收了！谢谢楼主！马后炮化工真好！

学会了 发表于 2010-4-16 16:46:47

精彩视频尽请收藏
网址1:d3eo7uqdxomcqx.cloudfront.net/02gjhttps://d3eo7uqdxomcqx.cloudfront.net/02gj
网址2：github.com/svappphttps://github.com/svappp

wunhua2004 发表于 2010-4-16 16:46:47


这东西我收了！谢谢楼主！马后炮化工真好！

jce 发表于 2010-4-16 16:46:47


论坛不能没有像楼主这样的人才啊！我会一直支持马后炮化工。

acrees 发表于 2010-4-16 16:46:47


论坛不能没有像楼主这样的人才啊！我会一直支持马后炮化工。

极地苍凉 发表于 2010-4-16 16:46:47

{:1110_550:}

王水 发表于 2024-8-28 10:11:54

{:1110_550:}
页: [1]
查看完整版本: 以前无聊的时候搞得个东东（二）