Fortran Coder

标题: 三次方程程序运行不了 [打印本页]

作者: muzili2008    时间: 2014-12-4 14:31
标题: 三次方程程序运行不了
大家好,我写了如下的代码,纠结三次方程。
[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
1、你这参数都不对应,怎么能编译;
作者: li913    时间: 2014-12-4 15:13
2、三次方程的根可能有复数,所以解X要用复数;
3、建议使用盛金公式,百度就有。你这个公式我没看懂。

作者: li913    时间: 2014-12-4 15:18
4、如果你只要实数解,修改  CALL CUBIC_ROOTS(IX,X)   就出结果了。
作者: muzili2008    时间: 2014-12-4 15:40
果然是这样,诶呀,大意了!!

多谢高手指教!!
作者: muzili2008    时间: 2014-12-4 16:41
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
对于复数赋值:
real a,b
complex c
c=cmplx(a,b)
作者: muzili2008    时间: 2014-12-5 15:53
li913 发表于 2014-12-5 15:29
对于复数赋值:
real a,b
complex c

原来如此 ,多谢了!!




欢迎光临 Fortran Coder (http://bbs.fcode.cn/) Powered by Discuz! X3.2