应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Windows编程 » 【求助】求助：vb编程中用牛顿迭代解三次方程为什么只得到一个根？【已完成】
52/1返回列表
查看: 2210  |  回复: 15
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
本帖产生 2 个 程序强帖 ，点击这里进行查看
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

liangzidou

银虫 (小有名气)

[交流] 【求助】求助：vb编程中用牛顿迭代解三次方程为什么只得到一个根？【已完成】已有5人参与

我在用vb编写中，遇到这样一个问题，解一个一元三次方程只得到一个根，实际上不只一个根。请各位高手帮忙找一下错误。
Private Sub Command1_Click()
Dim r, p, t, w, pc, tc, m, a, b, ca, cb, xl, xm, xn As Double
r = 8.314
p = 101.325
t = 301
pc = 3384
tc = 460.4
w = 0.227
m = 0.37464 + 1.5422 * w + 0.2699 * w ^ 2
tr = t / tc
a = 0.45723553 * r ^ 2 * tc ^ 2 * (1 + m * (1 - Sqr(tr))) ^ 2 / pc
b = 0.077796074 * r * tc / pc
Print "a="; a, "b="; b
ca = a * p / (r ^ 2 * t ^ 2)
cb = b * p / (r * t)
xl = cb - 1
xm = -3 * cb ^ 2 - 2 * cb + ca
xn = cb ^ 3 + cb ^ 2 - ca * cb
Print "A="; ca, "B="; cb
Print "xl="; xl, "xm="; xm
Print "xn="; xn
Dim f, z As Single
z = 0.01
Do
z = z1
f = z ^ 3 + xl * z ^ 2 + xm * z + xn
f1 = 3 * z ^ 2 + 2 * xl * z + xm
z1 = z - f / f1
Loop While Abs(z - z1) >= 0.0005
Print z1
End Sub

[ Last edited by liangzidou on 2010-5-15 at 22:22 ]
回复此楼

» 猜你喜欢
1楼2010-04-19 22:21:04
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

phychemlxd

金虫 (小有名气)
★ ★ ★ ★ ★ ★
小木虫(金币+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;
   }
赞一下(3人) 回复此楼
修.齐.治.平
11楼2010-04-21 18:00:29
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 16 个回答

liangzidou

银虫 (小有名气)
jjdg:把方程具体的值给出来 2010-04-20 12:52
方程就是f = z ^ 3 + xl * z ^ 2 + xm * z + xn=0
其中xl,xm,xn分别是方程的二次项、一次项系数及常数项。
下面这段程序就是我用迭代法就该方程的程序：

Dim f, z As Single
z = 0.01
Do
z = z1
f = z ^ 3 + xl * z ^ 2 + xm * z + xn
f1 = 3 * z ^ 2 + 2 * xl * z + xm
z1 = z - f / f1
Loop While Abs(z - z1) >= 0.0005
Print z1


最后得到的结果只有一个根。
赞一下(1人) 回复此楼
2楼2010-04-20 10:05:18
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

jjdg

版主 (知名作家)

优秀版主
★ ★
liangzidou(金币+2):谢谢了，我这个只能用vb程序解，因为下面的程序还要用到。 2010-04-20 12:37
wangen994(金币+2):辛苦了，jjdg 2010-04-20 14:13
不是还有其他解法吗？
我查到的其他解法是先判断delta与0的关系
①：当A=B=0时，方程有一个三重实根；
②：当Δ=B^2－4AC>0时，方程有一个实根和一对共轭虚根；
③：当Δ=B^2－4AC=0时，方程有三个实根，其中有一个两重根；
④：当Δ=B^2－4AC<0时，方程有三个不相等的实根。
赞一下(1人) 回复此楼
努力学习！以正当途径！获得需要的知识！
3楼2010-04-20 12:29:43
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

jjdg

版主 (知名作家)

优秀版主
我用盛金公式 写的
Private Sub Command1_Click()
Dim a, b, c, d, aa, bb, cc, delta, y1, y2, k, theta, t, x, x1, x2, x3, tt, tt1, xx2, xx3,s
'az^3+bz^2+cz+d=0，a、b、c、d来自一个例子（见最后）

a = 1
b = -70.5
c = 1533.54
d = -10082.44
aa = Round((b ^ 2 - 3 * a * c), 4)
bb = Round((b * c - 9 * a * d), 4)
cc = Round((c * c - 3 * b * d), 4)
delta = Round((bb * bb - 4 * aa * cc), 4)
If (aa = 0) And (bb = 0) Then
    x = -b / (3 * a)
    Label1.Caption = "X1=X2=X3=:" & x
Else
    Select Case delta
    Case Is > 0
        y1 = aa * b + 3 * a * (-bb + Sqr(delta)) / 2
        If y1 > 0 Then
            y1 = y1 ^ (1 / 3)
        Else
            y1 = y1 * (-1)
            y1 = y1 ^ (1 / 3) * (-1)
        End If
        y2 = aa * b + 3 * a * (-bb - Sqr(delta)) / 2
        If y2 > 0 Then
            y2 = y2 ^ (1 / 3)
        Else
            y2 = y2 * (-1)
            y2 = y2 ^ (1 / 3) * (-1)
        End If
        tt = y1 + y2
        tt1 = 3 ^ (1 / 3) * (y1 - y2) / (6 * a)
        s = (-2 * b + tt) / (6 * a)
        x1 = (-b - tt) / (3 * a)
        x2 = s & "+" & tt1 & "i"
        x3 = s & "-" & tt1 & "i"
        Label1.Caption = "X1=" & x1
        Label2.Caption = "X2=" & x2
        Label3.Caption = "X3=" & x3
    Case Is = 0
        k = bb / aa
        x1 = -b / a + k
        x = Str(-k / 2)
        Label1.Caption = "X1=" & x1
        Label2.Caption = "X2=" & x
        Label3.Caption = "X3=" & x
    Case Is < 0
        t = (2 * aa * b - 3 * a * bb) / (2 * aa ^ (3 / 2))
        theta = Atn(-t / Sqr(-t * t + 1)) + 2 * Atn(1)
        theta = theta / 3
        tt = aa ^ (1 / 2) * Cos(theta)
        tt1 = 3 ^ (1 / 2) * Sin(theta)
        x1 = (-b - 2 * tt) / (3 * a)
        xx2 = (-b + tt + tt1) / (3 * a)
        xx3 = (-b + tt - tt1) / (3 * a)
        Label1.Caption = "X1=" & x1
        Label2.Caption = "X2=" & xx2
        Label3.Caption = "X3=" & xx3        
    End Select
End If
End Sub

一建筑物的楼顶要建一个储水池，按施工的设计要求，这个储水池的长、宽、高之和为67.4dm，且宽=高，满储水量为9539.712(dm)＾3，立体对角线为1706.92dm，问：如何施工才能达到设计要求？
解：设取长、宽、高分别为X⑴、X⑵、X⑶，依题意：
X⑴＋X⑵＋X⑶=67.4；
X⑴X⑵X⑶=9539.712；
X⑴^2＋X⑵^2＋X⑶^2=1706.92。
解这个方程组，得一元三次方程
X＾3－67.4X＾2＋1417.92X－9539.712=0
a=1，b=－67.4，c=1417.92，d=－9539.712。
A=289；B=－9710.4；C=81567.36，
Δ=0。
根据盛金判别法，此方程有三个实根，其中两个相等。
应用盛金公式③求解。
K=—33.6。
把有关值代入盛金公式③，得：
X⑴=33.8(dm)；X⑵=X⑶=16.8(dm)。
经检验，结果正确。
因为宽=高，
所以，应取长为33.8dm；宽=高=16.8dm来进行施工。

[ Last edited by jjdg on 2010-4-20 at 14:20 ]
赞一下 回复此楼
努力学习！以正当途径！获得需要的知识！
4楼2010-04-20 14:18:56
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 16 个回答
普通表情 龙 兔 虎 猫 高级回复(可上传附件)
信息提示
关闭
请填处理意见
关闭