【求助】求助：vb编程中用牛顿迭代解三次方程为什么只得到一个根？【已完成】 - 程序语言 - Windows编程

版块导航: 正在加载中...

登录注册

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

24小时热门版块排行榜

返回列表

本帖产生 2 个程序强帖，点击这里进行查看

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

jjdg

版主 (知名作家)

程序强帖: 9
应助: 223 (大学生)
贵宾: 3.926
金币: 51003.4
红花: 238
帖子: 7227
在线: 474.2小时
虫号: 27000
注册: 2003-10-21
专业: 药物代谢与药物动力学
管辖: 程序语言

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
wangen994(金币+5):辛苦了 2010-04-20 17:38
liangzidou(金币+28):我之前没学过编程语言，现在学习热力学用到这些知识了，很是抓狂，多谢大侠帮忙，我试了试大侠给编的程序，解立方型方程很好用。太感谢了，辛苦了，呵呵。 2010-04-20 23:10
wangen994(金币+10, 程序强帖+1):辛苦了，呵呵 2010-04-21 08:30

或者：
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
'az^3+bz^2+cz+d=0
a = Val(Text1.Text)
b = Val(Text2.Text)
c = Val(Text3.Text)
d = Val(Text4.Text)

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
Private Sub Text1_Click()
Text1.SelStart = 0
Text1.SelLength = Len(Text1.Text)
End Sub
Private Sub Text2_Click()
Text2.SelStart = 0
Text2.SelLength = Len(Text2)
End Sub
Private Sub Text3_Click()
Text3.SelStart = 0
Text3.SelLength = Len(Text3)
End Sub
Private Sub Text4_Click()
Text4.SelStart = 0
Text4.SelLength = Len(Text4)
End Sub

[ Last edited by jjdg on 2010-4-20 at 14:40 ]

赞一下(2人)

回复此楼

努力学习！以正当途径！获得需要的知识！

5楼2010-04-20 14:34:06

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

智能机器人

Robot (super robot)

我们都爱小木虫

找到一些相关的精华帖子，希望有用哦~

点击这里搜索更多相关资源

科研从小木虫开始，人人为我，我为人人

phychemlxd

金虫 (小有名气)

程序强帖: 1
应助: 3 (幼儿园)
金币: 1219.1
红花: 8
帖子: 132
在线: 135.7小时
虫号: 230817
注册: 2006-03-26
专业: 化工热力学和基础数据

★ ★ ★ ★ ★ ★
小木虫(金币+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的回帖

相关版块跳转我要订阅楼主 liangzidou 的主题更新

返回列表