Aufgabe 6

PROGRAM B_09
 IMPLICIT NONE
 REAL, DIMENSION (3,3) :: A, p(0:3)
 REAL :: x
 INTEGER :: iostat, i, j, stellen
 CHARACTER*80 :: pol, px

 WRITE (*,'(A)') 'geben Sie die 3x3-Matrix zeilenweise ein: '
 DO i = 1, 3
  DO
   WRITE (*,'(A,I1,A)') 'Zeile ',i,':'
   READ (*,*,IOSTAT=iostat) (A(i,j),j=1,3)
   IF (iostat == 0) EXIT
  ENDDO
 ENDDO
!                                                     2      3
! characteristisches Polynom: det(A-xI) = p +p .x+p2.x + p .x
! Regel von Sarrus:                        0  1    2      3
! | a11-x  a12   a13  | a11-x  a12  |
! |  a21  a22-x  a23  |  a21  a22-x | = (a11-x).(a22-x).(a33-x) + a12.a23.a31 +
! |  a31   a32  a33-x |  a31   a32  |
!
! + a13.a21.a32 - a31.(a22-x).a13 - a32.a23.(a11-x) - (a33-x).a21.a12 =
!
! = (a11.a22-(a11+a22).x+x**2).(a33-x) + a12.a23.a31 + a13.a21.a32 -
! - a31.a22.a13 + a31.a13.x - a11.a32.a23 + a32.a23.x - a21.a12.a33 + a21.a12.x
!
! = a11.a22.a33 - (a11+a22).a33.x + a33.x**2 - a11.a22.x + (a11+a22).x**2 - x**3
! + a12.a23.a31 + a13.a21.a32 - a31.a22.a13 + a31.a13.x - a11.a32.a23 +
! + a32.a23.x - a21.a12.a33 + a21.a12.x =
!
! = -x**3 + (a11+a22).x**2 + a33.x**2 - 
! - (a11+a22).a33.x - a11.a22.x + a31.a13.x + a32.a23.x + a21.a12.x +
! + a11.a22.a33 + a12.a23.a31 + a13.a21.a32 - a31.a22.a13 - a11.a32.a23 -
! - a21.a12.a33
!
! = -x**3 + (a11+a22+a33).x**2 -
! - (a11.a33 + a22.a33 + a11.a22).x +
! + (a31.a13 + a32.a23 + a21.a12).x +
! + a11.a22.a33 + a12.a23.a31 + a13.a21.a32 -
! - (a31.a22.a13 + a11.a32.a23 + a21.a12.a33)

  p(0) = A(1,1)*A(2,2)*A(3,3) + A(1,2)*A(2,3)*A(3,1) + A(1,3)*A(2,1)*A(3,2) - &
        (A(3,1)*A(2,2)*A(1,3) + A(1,1)*A(3,2)*A(2,3) + A(2,1)*A(1,2)*A(3,3))
  p(1) = A(3,1)*A(1,3) + A(3,2)*A(2,3) + A(2,1)*A(1,2) - &
        (A(1,1)*A(3,3) + A(2,2)*A(3,3) + A(1,1)*A(2,2))
  p(2) = A(1,1)+A(2,2)+A(3,3)
  P(3) = -1
  WRITE (*,'(a)',ADVANCE='NO') 'Polynomkoeffizienten: ('
  DO j = 0, 3
   WRITE (pol,*) p(j)
   pol = ADJUSTL(pol)
   stellen = LEN_TRIM(pol)
   DO i = stellen, 1, -1
    IF (pol(i:i) == '0') THEN
     pol(i:i) = ' '
    ELSE
     EXIT
    ENDIF
   ENDDO
   WRITE (*,'(A)',ADVANCE='NO') TRIM(pol)
   IF (j /= 3) THEN
    WRITE (*,'(A)',ADVANCE='NO') ','
   ELSE
    WRITE (*,'(A)') ')'
   ENDIF
  ENDDO
 DO
  DO
   WRITE (*,'(A)',ADVANCE='NO') 'geben Sie x ein (Ende mit EOT): x = '
   READ (*,*,IOSTAT=iostat) x
   IF (iostat < 0) THEN
    STOP
   ELSE IF (iostat == 0) THEN
    EXIT
   ENDIF
  ENDDO
  WRITE (px,*) x**3*p(3) + x**2*p(2) + x*p(1) + p(0)
  WRITE (*,'(A)') 'ergibt p(x) = '//TRIM(ADJUSTL(px))
 END DO
END PROGRAM