Ask Question

Name:
Title:
Your Question:

Answer Question

Name:
Your Answer:
User Submitted Source Code!


Description:
  P1
Language: C/C++
Code:
      PROGRAM PRACTICA1
      IMPLICIT NONE

      INTEGER CONT, DIM
      PARAMETER(DIM=4)
      REAL X(DIM), Y(DIM), A0, A1

!     Inicializamos las variables

      OPEN(1, FILE='datos.dat')
      READ(1,*) (X(CONT), CONT=1, DIM)
      READ(1,*) (Y(CONT), CONT=1, DIM)
      CLOSE(1)
      
!     Llamamos a las subrutinas
      
      CALL MINCUAD(X, Y, DIM, A0, A1)
      CALL LAGRANGE(X, Y, DIM, 0.4, 0.8, 100)
      CALL NEWTON(X, Y, DIM, 0.4, 0.8, 100)
      
      END

!     ************************************************************************
!     SUBRUTINA MINCUAD:
!     Calcula los coeficientes de la recta de regresión por mínimos cuadrados
!     y los inserta en un fichero (mincuad.dat)
!     RECIBE: Los vectores de abscisas (X) y ordenadas (Y), su dimensión (DIM)
!     ************************************************************************
      
      SUBROUTINE MINCUAD(X, Y, DIM, A0, A1)
      INTEGER DIM
      REAL X(DIM), Y(DIM), A0, A1

      INTEGER CONT
      DOUBLEPRECISION SUMAX, SUMAX2, SUMAXY, SUMAY
      
!     Calculamos los sumatorios

      SUMAX=0
      SUMAX2=0
      SUMAXY=0
      SUMAY=0
      DO CONT=1, DIM
         SUMAX=SUMAX+X(CONT)
         SUMAX2=SUMAX2+X(CONT)**2
         SUMAXY=SUMAXY+X(CONT)*Y(CONT)
         SUMAY=SUMAY+Y(CONT)
      ENDDO

!     Calculamos los coeficientes
      
      A0=(SUMAY*SUMAX2-SUMAX*SUMAXY)/(DIM*SUMAX2-SUMAX**2)
      A1=(DIM*SUMAXY-SUMAY*SUMAX)/(DIM*SUMAX2-SUMAX**2)

      OPEN(1, FILE='mincuad.dat')
      WRITE(1,*) 'y=a0+a1x'
      WRITE(1,*) 'a0=', A0
      WRITE(1,*) 'a1=', A1
      CLOSE(1)

      RETURN
      END

!     ************************************************************************
!     SUBRUTINA LAGRANGE:
!     Calcula el polinomio interpolador de Lagrange y evalúa una serie de
!     puntos en él, escribiéndolos en un fichero (Lagrange.dat)
!     RECIBE: los vectores de abscisas (X) y de ordenadas (Y), su dimensión
!     (DIM), el intervalo (INI, FIN) y número de puntos a interpolar
!     (NUMPUNTOS)
!     ************************************************************************

      SUBROUTINE LAGRANGE(X, Y, DIM, INI, FIN, NUMPUNTOS)
      INTEGER DIM, NUMPUNTOS
      REAL X(DIM), Y(DIM), INI, FIN

      INTEGER I, J
      REAL STEP, VAR, FINNUEVO
      DOUBLEPRECISION SUMA, PROD

!     Calculamos el incremento de la variable VAR
      
      STEP=(FIN-INI)/(1.0*NUMPUNTOS)
      FINNUEVO=FIN+STEP
      
      
      OPEN(1, FILE='Lagrange.dat')

!     Obtenemos el valor numérico del polinomio en cada punto
      
      VAR=INI
      DO WHILE(VAR.LE.FINNUEVO)
         SUMA=0        
         DO I=1, DIM
            PROD=Y(I)
            DO J=1, I-1
               PROD=PROD*(VAR-X(J))/(X(I)-X(J))
            ENDDO
            DO J=I+1, DIM
               PROD=PROD*(VAR-X(J))/(X(I)-X(J))
            ENDDO
            SUMA=SUMA+PROD
         ENDDO
         WRITE(1,*) 'f(', VAR, ')=', SUMA
         VAR=VAR+STEP
      ENDDO
      
      CLOSE(1)

      RETURN
      END

!     ************************************************************************
!     SUBRUTINA NEWTON:
!     Calcula el polinomio interpolador de Newton y evalúa una serie de
!     puntos en él, escribiéndolos en un fichero (Newton.dat)
!     RECIBE: los vectores de abscisas (X) y de ordenadas (Y), su dimensión
!     (DIM), el intervalo (INI, FIN) y número de puntos a interpolar
!     (NUMPUNTOS)
!     ************************************************************************

      SUBROUTINE NEWTON(X, Y, DIM, INI, FIN, NUMPUNTOS)
      INTEGER DIM, NUMPUNTOS
      REAL X(DIM), Y(DIM), INI, FIN

      INTEGER I, J
      REAL STEP, VAR, FINNUEVO
      DOUBLEPRECISION DIFDIV(DIM, DIM), SUMA, PROD

!     Calculamos las diferencias divididas
      
      DO I=1, DIM
         DIFDIV(I, 1)=Y(I)
         DO J=2, I
            DIFDIV(I, J)=(DIFDIV(I, J-1)-DIFDIV(I-1, J-1))/(X(I)
     &           -X(I-J+1))
         ENDDO
      ENDDO

      STEP=(FIN-INI)/(1.0*NUMPUNTOS)
      FINNUEVO=FIN+STEP

      OPEN(1, FILE='Newton.dat')
      
!     Obtenemos el valor numérico del polinomio en cada punto
      
      VAR=INI
      DO WHILE(VAR.LE.FINNUEVO)
         SUMA=0
         DO I=1, DIM
            PROD=DIFDIV(I, I)
            DO J=1, I-1
               PROD=PROD*(VAR-X(J))
            ENDDO
            SUMA=SUMA+PROD
         ENDDO
         WRITE(1,*) 'f(', VAR, ')=', SUMA
         VAR=VAR+STEP
      ENDDO
      
      CLOSE(1)

      RETURN
      END
        
Comments: