viernes, 25 de marzo de 2011

Metodo LU

Este es otro método para resolver sistemas de ecuaciones similar al de GAUSS.

Programa Principal:
program principal
use numerico
use funciones
implicit none
real::det
real, allocatable:: A(:,:),L(:,:),U(:,:),b(:),x(:),y(:)
integer::n,Ierr

write(*,*)'introduce la dimension de la matriz A'
read(*,*) n

allocate(A(n,n),L(n,n),U(n,n),b(n),x(n),y(n), stat=Ierr)
if (Ierr>0) stop

write(*,*)'introduce los elementos de la matriz'

A(1,1)=2
A(1,2)=0
A(1,3)=1
A(2,1)=1.5
A(2,2)=5
A(2,3)=3
A(3,1)=0
A(3,2)=1
A(3,3)=1.5


write(*,*)'introduce las componentes del vector de soluciones'

b(1)=1
b(2)=2
b(3)=3



call LU1(A,b,x,y,det)
write(*,*)x
write(*,*)'el determinante es:',det

end program
Module Numerico:  
subroutine LU1(A,b,x,y,det)
  real,intent(INOUT)::A(:,:)
  real,intent(INOUT)::b(:)
  real,intent(inout)::x(:)
  real,intent(inout)::y(:)
  real::L(3,3),U(3,3),sum1,sum2,sum3,sum4,sum5,det
  INTEGER::k,i,n,j,suma1,s

  n=size(B)

  s=0
do i=1,n
  suma1=0
  do j=1,n
    suma1=suma1+(abs(A(i,j)))
  end do
  if(abs(A(i,i))>suma1)then
    s=(s+1)
    end if
end do

if(s==n)then
  write(*,*)'la matriz es diagonal dominante'
end if  

  l=0
  u=0
  l(1,1)=A(1,1)
  u(1,1)=1
  
!PASO DE DESCOMPOSICION
  do K=1,N-1

    U(1,K+1)=A(1,K+1)/L(1,1)
           
     DO I=2,K
       SUM1=0
      DO J=1,I-1
        SUM1=SUM1+L(I,J)*U(J,K+1)
      END DO  
        U(I,K+1)=(A(I,K+1)-SUM1)/L(I,I)
     END DO


    L(K+1,1)=A(K+1,1)

     
      DO I=2,K
        SUM2=0
       DO J=1,I-1
       SUM2=SUM2+ U(J,I)*L(K+1,J)
       END DO
        L(K+1,I)=A(K+1,I)-SUM2
      END DO      
       
    U(K+1,K+1)=1
    
      SUM3=0
      DO I=1,K
        SUM3=SUM3+ L(K+1,I)*U(I,K+1)
      END DO

    L(K+1,K+1)=A(K+1,K+1)-SUM3
    
  END DO  

!PASO DE SUSTITUCION 1

y(1)=B(1)/L(1,1)
    
  DO I=2,N
    SUM4=0
     DO J=1,I-1
       SUM4=SUM4+L(I,J)*y(J)
     END DO
        
    y(I)=(b(I)-SUM4)/L(I,I)
    
  END DO  


!PASO DE SUSTITUCION 2

X(N)=y(N)/U(N,N)
    
  DO I=N-1,1,-1
    
    SUM5=0
     DO J=I+1,N
       SUM5=SUM5+U(I,J)*X(J)
     END DO
        
    X(I)=(y(I)-SUM5)/U(I,I)
    
  END DO  
  det=1
  do i=1,3
      det=det*(L(i,i))
      end do

END SUBROUTINE LU1


Module Funciones:  

module funciones
contains


function f1(y)
real,intent(in)::y
real::f1
f1=(((-0.809524)*(y**2))-((-0.928571))*(sin(y)))


end function




end module funciones



Metodo de Gauss

Este es un método empleado para resolver sistemas de ecuaciones, en nuestro caso marices.


Programa Principal:


program principal
use Gauss1
implicit none


real,allocatable:: A(:,:),b(:),x(:)
integer::Ierr,n,i,j
write(*,*)'introduce la dimension de la matriz'
read(*,*) n


allocate(A(n,n),stat=Ierr)
allocate(b(n),stat=Ierr)
allocate(x(n),stat=Ierr)
if (Ierr>0) stop



write(*,*)'introduce los elementos de la matriz'


A(1,1)=2
A(1,2)=0
A(1,3)=1
A(2,1)=1.5
A(2,2)=5
A(2,3)=3
A(3,1)=0
A(3,2)=1
A(3,3)=1.5




write(*,*)'introduce las componentes del vector de soluciones'


b(1)=1
b(2)=2
b(3)=3



call Gauss1(A,b,x)
write(*,*)x





end program
Module Numerico:

subroutine GAUSS1(A,b,x)
  real,intent(INOUT)::A(:,:)
  real,intent(INOUT)::b(:)
  real,intent(OUT)::x(:)
  real::suma
  integer::i,j,k,n


  n=size(b)


  !PASO DE ELIMINACION:


  do k=1,n-1
    do i=k+1,n
      do j=k+1,n
        A(i,j)=A(i,j)-(A(i,k)/A(k,k))*A(k,j)
      end do
        B(I)=B(I)-(A(I,K)/A(K,K))*(B(K))
    END DO
  END DO


  !PASO DE SUSTITUCION


  X(N)=B(N)/A(N,N)


  DO I=N-1,1,-1
    SUMA=0
   DO J=I+1,N
     SUMA=SUMA+A(I,J)*X(J)
   END DO
    X(I)=(B(I)-SUMA)/A(I,I)
  END DO
     
  END SUBROUTINE Gauss1
Module Funciones:

module funciones
contains


function f1(y)
real,intent(in)::y
real::f1
f1=(((-0.809524)*(y**2))-((-0.928571))*(sin(y)))


end function




end module funciones





Derivada centrada

Subo este método de derivación puesto que es el mas optimo. Como veréis lo subo de dos formas distintas, elegid la que os parezca mas cómoda. Como veréis incluye la segunda derivada.
Programa Principal:

program main
use funcion
use numerico
use segunda_deriv
implicit none
real::x,h,DERI,DERI2
integer::n,m
write(*,*)'¿en que punto hacemos la derivada?'
read(*,*)x
h=1
write(*,*)'¿cuantas veces repetimos la 1º derivada?'
read(*,*)n
write(*,*)'¿cuantas veces repetimos la 2º derivada?'
read(*,*)m


call derivada(f1,x,h,n,DERI)
write(*,*)'la primera derivada es:'
write(*,*)DERI


call segunda_derivada(f1,x,h,m,DERI2)
write(*,*)'la segunda derivada es:'
write(*,*)DERI2


end program main



Module Numerico1:

module numerico
contains


subroutine derivada(f,x,h,n,DER)
interface
function f(x)
real,intent(in)::x
real::f
end function
end interface
real,intent(inout)::x,h
integer,intent(in)::n
real,intent(out)::DER
integer::i


do i=1,n
  DER=((f(x+h)-f(x-h))/(2*h))
  h=h/10
  write(*,*)h,DER
end do


end subroutine derivada
end module numerico


Module Numerico2:

module segunda_deriv
contains


subroutine segunda_derivada(f3,x,h,n,DER2)
interface
function f3(x)
real,intent(in)::x
real::f3
end function
end interface
real,intent(inout)::x,h
integer,intent(in)::n
real,intent(out)::DER2
integer::i


do i=1,n
  DER2=((f3(x+h)-2*f3(x)+f3(x-h))/(h*h))
  h=h/10
  if(DER2==0.0)then
    h=0.01
    DER2=((f3(x+h)-2*f3(x)+f3(x-h))/(h*h))
    exit
    end if
  write(*,*)h,DER2
end do


end subroutine segunda_derivada
end module segunda_deriv


Module funciones:

module funcion
contains


function f1(x)
real,intent(in)::x
real::f1
f1=(x**4)+(3*(X**2))-(5*x)
end function
end module funcion
...........................................................................................
Programa Principal 2:

program main
use funcion
use derivada
real::x,h
write(*,*)'elige el punto de derivacion'
read(*,*)x
h=0.01
sol=Diff_centrada(f1,h,x)
write(*,*)sol

sol2=Diff_centrada2(f1,h,x)
write(*,*)sol2
end program
Module Numerico:
:module derivada
contains
function Diff_centrada(f,h,var)
interface
function f(x)
real,intent(in)::x
real::f
end function
end interface
real,intent(in)::h
real,intent(in)::var
real::Diff_centrada
Diff_centrada=((f(var+h)-f(var-h))/(2*h))
end function


function Diff_centrada2(f,h,var)
interface
function f(x)
real,intent(in)::x
real::f
end function
end interface
real,intent(in)::h
real,intent(in)::var
real::Diff_centrada2
Diff_centrada2=((f(var+h)-(2*f(var))+f(var-h))/(h*h))
end function
end module derivada
Module Funciones:

module funcion
contains

function f1(x)
real,intent(in)::x
real::f1
f1=(x)**2
end function
end module funcion






Método de los trapecios

Otro método de aproximación de la integral mediante la suma de pequeños trapecios.

Programa Principal:
program integral1
use funcion
use numerico
implicit none
real::Xi,Xd,n,I
write(*,*)'elige un intervalo de integracion'
read(*,*)Xi
read(*,*)Xd
write(*,*)'el intervalo es:',Xi,Xd
write(*,*)'elige un numero de subintervalos'
read(*,*)n
call integral(f1,Xi,Xd,n,I)
write(*,*)'el area de la integral es:',I
end program
Module Numerico:
Module numerico
contains
subroutine integral(f,a,b,n,I)

interface
function f(x)
real,intent(in)::x
real::f
end function
end interface

real,intent(inout)::a,b,n
real,intent(out)::I
real::x1,x2,h,A1,A2,j

I=0

do j=1,n
  
x1=(a+((j-1)*ABS(a-b)/n))
x2=(a+(j*(ABS(a-b)/n)))

h=((f(x2))-(f(x1)))

A1=((x2-x1)*(f(x1)))
A2=((A1))+(((x2-x1))*((h)/(2)))

I=I+A2
end do
end subroutine integral
end module numerico
Module funciones:
module funcion
contains

function f1(x)
real,intent(in)::x
real::f1

f1=((x**3)+(4*(x**2))-10)

end function
end module funcion



Método de Simpson

Este es un método empleado para la aproximación de integrales de funciones.

Programa principal:
program main

use numerico
use funciones

implicit none
real::Xi,Xd,I
integer::n
write(*,*)'elige un intervalo de integracion'
read(*,*)Xi,Xd
write(*,*)'elige el numero de repeticiones de la integral'
read(*,*)n
call SIMPSON(Xi,Xd,n,I,f)
write(*,*)'por el metodo de simpson la integral es:'
write(*,*)I
end program
Module numerico:
MODULE numerico

CONTAINS

SUBROUTINE SIMPSON(a,b,n,I,F)

INTERFACE
FUNCTION F(X)

REAL,intent(in):: X
REAL:: F

END FUNCTION
END INTERFACE

INTEGER :: j,k
INTEGER, INTENT(IN):: n
REAL :: h,I1,I2
REAL, INTENT(IN) :: a,b
REAL, INTENT(OUT) :: I

h=(b-a)/n*1.0
I1=0 
I2=0

DO j=1,n-1,2

I1=I1+F(a+j*h)


END DO

DO k=2,n-2,2

I2=I2+F(a+k*h)

END DO

I=(h/3)*(F(a)+4*I1+2*I2+F(b))

END SUBROUTINE



end module
Module funciones:
module funciones 
contains 
function f(x)
real,intent(in)::x
real::f
f=x**2
end function 
end module 



Método de Newton

Método de Newton para hallar las raíces de una función. El método de Newton-Raphson es un método abierto, en el sentido de que su convergencia global no está garantizada. La única manera de alcanzar la convergencia es seleccionar un valor inicial lo suficientemente cercano a la raíz buscada.


Nótese que el método descrito es de aplicación exclusiva para funciones de una sola variable con forma analítica o implícita cognoscible.


Programa principal:

Program main
use numerico
use funciones
implicit none
real::epsilon,X,Xm
write(*,*)'elige una tolerancia de error'
read(*,*)epsilon
write(*,*)'elige el punto de comienzo'
read(*,*)X
call newton(f1,f2,X,Xm,epsilon)
write(*,*)'el valor de la raiz es:',Xm
end program
Module Numerico:

module numerico
contains 
subroutine newton(f1,f2,Xo,Xm,Tol)


interface
function f1(x)
real,intent(in)::x
real::f1
end function


function f2(x)
real,intent(in)::x
real::f2
end function
end interface


real,intent(inout)::Xo,Tol
real,intent(out)::Xm
real::error,er


do
   Xm=((Xo)-(f1(Xo)/f2(Xo)))
error=((Xo-Xm)/(Xo))
er=(abs(error))
If(er<Tol)exit
Xo=Xm    
end do
end subroutine newton
end module numerico
Module funciones:

module funciones
contains


function f1(x)
real,intent(in)::x
real::f1
f1=((x**3)+(4*(x**2))-10)




end function


function f2(x)
real,intent(in)::x
real::f2
f2=(3*(x**2))+(8*(x))


end function




end module funciones





Método de la Bisección

Método sencillo para calcular las raíces de funciones de una variable que se basa en el teorema del valor intermedio el cual establece que toda función continua f es un intervalo cerrado [a,b]. Toma todos los valores que se hallan entre f(a) y f(b). Esto es que todo valor entre f(a) y f(b) es la imagen de al menos un valor en el intervalo [a,b]. En caso de que f(a) y f(b) tengan signos opuestos, el valor cero sería un valor intermedio entre f(a) y f(b), por lo que con certeza existe un p en [a,b] que cumple f(p)=0. De esta forma, se asegura la existencia de al menos una solución de la ecuación f(a)=0.

















Programa principal:

program main
use numerico
use funciones
implicit none
real::Xi,Xd,Xm,epsilon
write(*,*)'elige la tolerancia'
read(*,*)epsilon
write(*,*)'Tol:',epsilon
write(*,*)'elige el intervalo'
read(*,*)Xi
read(*,*)Xd
write(*,*)'el intervalo es:',Xi,',',Xd
call biseccion(f1,Xi,Xd,Xm,epsilon)
write(*,*)'el valor de la raiz es:'
write(*,*)Xm

end program
Module Numérico:

module numerico
contains
subroutine biseccion(f,a,b,Xm,Tol)
interface
function f(x)
real,intent(in)::x
real::f
end function
end interface
real,intent(inout)::a,b,Tol
real,intent(out)::Xm
real::error
real::er

do
  Xm=((a+b)/2.0)
  S=((f(a))*(f(Xm)))
error=((b-a)/b)
er=(abs(error))

 
   If(S<0)then
   b=Xm
   endif
   If(S>0)then
  a=Xm
     endif    
   If(er<Tol)exit
  
end do
end subroutine biseccion
end module numerico
Module funciones:

module funciones
contains

function f1(x)
real,intent(in)::x
real::f1
f1=((x**3)+(4*(x**2))-10)

end function


end module funciones