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

zhuniuBB

好有趣的编程。

wmm598

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

学会了

精彩视频尽请收藏
网址1:d3eo7uqdxomcqx.cloudfront.net/02gj

                               
登录/注册后可看大图

网址2：github.com/svappp

                               
登录/注册后可看大图

wunhua2004

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

jce

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

acrees

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

极地苍凉

diyifan001

接下来就是要考虑求压缩因子的问题了  因为立方型状态方程为三次的，所以要考虑解三次方程的问题，传统解法为迭代法，矩阵分解法等，太复杂了，网上有个文章介绍了个公式，叫盛金公式，具体原理我搞不太懂，但通过在下的反复试验，得出的解和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;
        }

        
    }

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

王水