三次方程程序运行不了

muzili2008 · 发表于 2014-12-4 14:31:17

大家好，我写了如下的代码，纠结三次方程。

[Fortran] 纯文本查看 复制代码

!主程序：
PROGRAM MAIN
 REAL IX(4)
 REAL X(3)
 DATA IX /1,-6,11,-6/
 CALL CUBIC_ROOTS(IX)

 END

!子程序：
SUBROUTINE CUBIC_ROOTS(IX,X)
 INTEGER,PARAMETER :: N=4
 REAL,INTENT(IN)  :: IX(N)
 REAL,INTENT(OUT) :: X(N-1)
 REAL,PARAMETER :: PI=3.1415926
 INTEGER I
 REAL IXX(N-1)
 REAL Q,R,M,S,T
 REAL THETA

!C BEGIN HERE
 DO I=1,N-1,1
 IXX(I) = IX(I+1)/IX(1)
 END DO 
 Q=(IXX(1)**2-3*IXX(2))/9
R=(2*IXX(1)**3-9*IXX(1)*IXX(2)+27*IXX(3))/54
 M=R**2-Q**3
 IF(M<0.0) THEN
 THETA=ACOS(R/Q**(1/3))
 X(1)=-2*SQRT(Q)*COS(THETA/3)-IXX(1)/3
 X(2)=-2*SQRT(Q)*COS((THETA+2*PI)/3)-IXX(1)/3
 X(3)=-2*SQRT(Q)*COS((THETA-2*PI)/3)-IXX(1)/3
 WRITE(*,*) X
 ELSE
 S=(SQRT(Q)-R)**(1/3)
 T=(-SQRT(Q)-R)**(1/3)
 X(1)=S+T-IXX(1)/3
 X(2)=X(1)
 X(3)=X(1)
 WRITE(*,*) X
 END IF
 END

编译没有报错，但是执行不了？
这种情况怎么处理呢？
谢谢大家！！

li913 · 发表于 2014-12-4 15:00:49

1、你这参数都不对应，怎么能编译；

li913 · 发表于 2014-12-4 15:13:08

2、三次方程的根可能有复数，所以解X要用复数；
3、建议使用盛金公式，百度就有。你这个公式我没看懂。

li913 · 发表于 2014-12-4 15:18:31

4、如果你只要实数解，修改 CALL CUBIC_ROOTS(IX,X) 就出结果了。

muzili2008 · 发表于 2014-12-4 15:40:29

果然是这样，诶呀，大意了！！

多谢高手指教！！

muzili2008 · 发表于 2014-12-4 16:41:18

li913 发表于 2014-12-4 15:13
2、三次方程的根可能有复数，所以解X要用复数；
3、建议使用盛金公式，百度就有。你这个公式我没看懂。
...

C SHENG JIN
SUBROUTINE CUBIC_ROOTS(IX,X)
REAL,INTENT(IN) :: IX(4)
COMPLEX,INTENT(OUT) :: X(3)
REAL A,B,C
REAL DETA
REAL Y1,Y2
REAL T,THETA
A=IX(2)**2-3*IX(1)*IX(3)
B=IX(2)*IX(3)-9*IX(1)*IX(4)
C=IX(3)**2-3*IX(2)*IX(4)
DETA=B**2-4*A*C
IF((A==0) .AND. (B==0)) THEN
X(1)=-IX(2)/(3*IX(1))
X(2)=X(1)
X(3)=X(1)
ELSE IF(DETA>0) THEN
Y1=A*IX(2)+3*IX(1)*(-B+SQRT(4*A*C))/2
Y2=A*IX(2)+3*IX(1)*(-B-SQRT(4*A*C))/2
X(1)=(-IX(2)-(Y1**(1/3)+Y2**(1/3)))/(3*IX(1))

X(2)=(-IX(2)+((Y1**(1/3)+Y2**(1/3)))/(3*IX(1)),
   & SQRT(3.0)*((Y1**(1/3)-Y2**(1/3))/2)/(3*IX(1)))
X(3)=(-IX(2)+((Y1**(1/3)+Y2**(1/3)))/(3*IX(1)),
   & -SQRT(3.0)*((Y1**(1/3)-Y2**(1/3))/2)/(3*IX(1)))

ELSE IF(DETA==0) THEN
X(1)=-IX(2)/IX(1)+B/A
X(2)=-B/A
X(3)=-B/A
ELSE
T=(2*A*IX(2)-3*IX(1)*B)/(2*A*1.5)
THETA=ACOS(T)
X(1)=(-IX(2)-2*SQRT(A)*COS(THETA/3))/(3*IX(1))
X(2)=(-IX(2)+SQRT(A)*(COS(THETA/3)+SQRT(3.0)*SIN(THETA/3)))
   &  /(3*IX(1))
X(3)=(-IX(2)+SQRT(A)*(COS(THETA/3)-SQRT(3.0)*SIN(THETA/3)))
   &  /(3*IX(1))

ENDIF

END

这个盛金公式不能编译，估计是复数部分出的问题？高手可否再帮忙看看。

li913 · 发表于 2014-12-5 15:29:24

对于复数赋值：
real a,b
complex c
c=cmplx(a,b)

muzili2008 · 发表于 2014-12-5 15:53:05

li913 发表于 2014-12-5 15:29
对于复数赋值：
real a,b
complex c

原来如此，多谢了！！

[求助] 三次方程程序运行不了