Tuesday, November 16, 2010

Fortran program to calculate an integral using "crude" Monte Carlo

The following program calculates the integral I

using "crude" Monte Carlo. In this case, n=3, xia=0, xib=1 and f(x)=exp(a·x).


Program intmc

double precision fsum,fsum2,a,xa,xb,V,x,f,I,IA,varf
integer M,n,j,k
dimension a(3),xa(3),xb(3),x(3)
c n: dimension of the array x
c M: number of random numbers we use
c V: volum of the hipervolume
n=3
M=1000000
V=1
do j=1,n
a(j)=rand()
xa(j)=0
xb(j)=1
V=V*(xb(j)-xa(j))
enddo
c let's find the value of the integral that we get using monte-carlo
fsum=0
c fsum2 will help us to find the error
fsum2=0
do k=1,M
do j=1,n
x(j)=rand()
enddo
call func(a,x,f,n)
fsum=fsum+f
fsum2=fsum2+f**2
enddo
c I is the value of the integral we wanted to find
c IA is the value of the integral calculated analiticaly
I=V*fsum/M
varf=(fsum2/M-((fsum/M)**2))
sigma=sqrt(varf/M)
IA=1
do j=1,n
IA=IA*((exp(a(j))-1)/a(j))
enddo
write(*,*) I,sigma*3,IA
stop
end


SUBROUTINE func(a,x,f,n)
double precision a,x,f,ax
integer n
dimension a(n),x(n)
ax=0
do j=1,n
ax=ax+a(j)*x(j)
enddo
f=exp(ax)
end

Fortran program to calculate an integral using the method of trapezes

The following program calculates the integral I

using the trapezes method. In this case, n=3, xia=0, xib=1 and f(x)=exp(a·x).


Program inttr
double precision isum,x,xa,xb,h,f,I,a,IA,error,e
integer j,p,k1,k2,k3,n,l
dimension x(3),xa(3),xb(3),h(3),a(3)
n=3
p=300
do j=1,n
a(j)=rand()
xa(j)=0
xb(j)=1
h(j)=(xb(j)-xa(j))/p
enddo
isum=0
error=0
c we have to use n do's one inside the other. In this case, n=3.
c to calculate the error we need to compute the second partial
c derivate of F
do k1=1,p
x(1)=xa(1)+(k1-1/2)*h(1)
do k2=1,p
x(2)=xa(2)+(k2-1/2)*h(2)
do k3=1,p
x(3)=xa(3)+(k3-1/2)*h(3)
call func(a,x,f,n)
isum=isum+f
do l=1,n
error=error+f*(a(l)**2)
enddo
enddo
enddo
enddo
I=isum*(h(1)*h(2)*h(3))
e=(p**(-2/n))*error/(24*isum)
IA=1
do j=1,n
IA=IA*((exp(a(j))-1)/a(j))
enddo
write(*,*) I,e,IA
stop
end



SUBROUTINE func(a,x,f,n)
double precision a,x,f,ax
integer n,j
dimension a(n),x(n)
ax=0
do j=1,n
ax=ax+a(j)*x(j)
enddo
f=exp(ax)
end

Metropolis Algorithm

A very interesting link where we can understand easily the Metropolis Algorithm:

http://www.childrens-mercy.org/stats/weblog2007/MetropolisAlgorithm.asp

Saturday, November 6, 2010

Calculate PI with a simple fortran program

In this link (Fer problemes), you can find a montecarlo problem that allows us to find the number PI.

Here, you have one code that could be implemented to simulate this program. It is done in fortran.

PROGRAM calculpi

integer b,d,R
real acirc, aquad
double precision div
b=12345
d=76545
R=1
acirc=0
aquad=0
do i=1,1000000
a=rand(b)
c=rand(d)
b=a
d=c
if ((a**2+c**2)
acirc=acirc+1
endif
if ((a
aquad=aquad+1
endif
div=(acirc/aquad)*4

if (mod(i,10)==0)then
write (*,*) acirc, aquad, div
endif
enddo
STOP
END

Friday, November 5, 2010

3D plots in gnuplot

The comand is:

gnuplot> splot "positions.dat"

Wednesday, November 3, 2010

Problem in fortran using lapack

Consider the linear system of equations Ax=b with the block structure:

We know that the LU decomposition of A is:

Make a program in fortran to compute the LU decomposition and to find the solution.

_________________________________________________________________



Program ludecomposition

integer INFO1,INFO2, i, j, N, N1,N2,N3
integer, ALLOCATABLE, dimension (:) :: IPIVZ1,IPIVZ2
double precision, ALLOCATABLE, DIMENSION(:,:) :: A11,A13,A22,
c A23,A31,A32,A33,S,Z3,Z4,A11C,A22C
double precision, ALLOCATABLE, DIMENSION(:) ::B1,B2,B3
c we let the user to enter the dimensions N1,N2 and N3.
open (unit=1, file="solution.dat")
write(*,*) 'Input the 3 values n1,n2,n3.After each press enter.'
read (*,*) N1,N2,N3
c The dimension of A is N. Let's compute N.
N=N1+N2+N3
C A is a square matrix. That's why we have to compute M=N.
ALLOCATE(A11(N1,N1),A13(N1,N3),A22(N2,N2),A23(N2,N3),
c A31(N3,N1),A32(N3,N2),A33(N3,N3),S(N3,N3),Z3(N3,N1),
c Z4(N3,N2),A11C(N1,N1),A22C(N2,N2))
ALLOCATE(B1(N1),B2(N2),B3(N3),IPIVZ1(N1),IPIVZ2(N2))
C Let's full the matrixes AII with random numbers.
do i=1,N1
do j=1,N1
A11(i,j)=rand()
enddo
do j=1,N3
A13(i,j)=rand()
enddo
enddo
do i=1,N2
do j=1,N2
A22(i,j)=rand()
enddo
do j=1,N3
A23(i,j)=rand()
enddo
enddo
do i=1,N3
do j=1,N1
A31(i,j)=rand()
enddo
do j=1,N2
A32(i,j)=rand()
enddo
do j=1,N3
A33(i,j)=rand()
enddo
enddo
C Now the matrixs AII are full of random numbers. In order to see if
C it works correcly, we can full the matrix b with the elements of
C each row added.
do i=1,N1
B1(i)=0
do j=1,N1
B1(i)=B1(i)+A11(i,j)
enddo
do j=1,N3
B1(i)=B1(i)+A13(i,j)
enddo
enddo
do i=1,N2
B2(i)=0
do j=1,N2
B2(i)=B2(i)+A22(i,j)
enddo
do j=1,N3
B2(i)=B2(i)+A23(i,j)
enddo
enddo
do i=1,N3
B3(i)=0
do j=1,N1
B3(i)=B3(i)+A31(i,j)
enddo
do j=1,N2
B3(i)=B3(i)+A32(i,j)
enddo
do j=1,N3
B3(i)=B3(i)+A33(i,j)
enddo
enddo
C now, we call our subroutine in order to compute the LU decomposition
C of the matrix A composed by our AII matrixes
call ludec(A11,A13,A22,A23,A31,A32,A33,
C N1,N2,N3,S,Z3,Z4,A11C,A22C,IPIVZ1,IPIVZ2)
C finally, we call our subrouine in order to get the values of x
call solve(A11,A13,A22,A23,A31,A32,A33,
C N1,N2,N3,S,Z3,Z4,B1,B2,B3,A11C,A22C,IPIVZ1,IPIVZ2)
C let's write the result of the computation. If we write B, we are
C writting the solution of the system.
write(1,*)'The solution of the x-array you wanted to know is'
do i=1,N1
write(1,*)B1(i)
enddo
do i=1,N2
write(1,*)B2(i)
enddo
do i=1,N3
write(1,*)B3(i)
enddo
write(*,*)'The solution can be found in the file solution.dat.'
C in our example, x-array should be a vector of ones.
stop
end

C==================================================
subroutine ludec(A11,A13,A22,A23,A31,A32,A33,
C N1,N2,N3,S,Z3,Z4,A11C,A22C,IPIVZ1,IPIVZ2)
integer INFO,i,j,k
CHARACTER(1) TRANSN,TRANST
integer IPIVZ1(N1),IPIVZ2(N2)
double precision A11,A13,A22,A23,A31,A32,A33,
C A31Z1,A32Z2,S,A31T,A32T,Z3,Z4,Z1,Z2,A11C,A22C
dimension A11(N1,N1),A13(N1,N3),A22(N2,N2),A23(N2,N3),
C A31(N3,N1),A32(N3,N2),A33(N3,N3),A(N1+N2+N3,N1+N2+N3),
C A31Z1(N3,N3),A32Z2(N3,N3),S(N3,N3),A31T(N1,N3),Z3(N3,N1),
C A32T(N2,N3),Z4(N3,N2),Z1(N1,N3),Z2(N2,N3),
C A11C(N1,N1),A22C(N2,N2)
TRANSN='N'
TRANST='T'
N=N1+N2+N3
C let's make a copy of A11,A22,A13 and A23
do i=1,N1
do j=1,N1
A11C(i,j)=A11(i,j)
enddo
do j=1,N3
Z1(i,j)=A13(i,j)
enddo
enddo
do i=1,N2
do j=1,N2
A22C(i,j)=A22(i,j)
enddo
do j=1,N3
Z2(i,j)=A23(i,j)
enddo
enddo
c let's find the value of Z1=(A11-1)*A13. We know that A11*Z1=A13
c let's get the LU decomposition of A11
call DGETRF( N1, N1, A11C, N1, IPIVZ1, INFO)
c let's solve the system A11*Z1=A13 to find Z1.
call DGETRS(TRANSN,N1,N3,A11C,N1,IPIVZ1,Z1,N1,INFO)
c=====
c let's find the value of Z2=(A22-1)*A23. We know that A22*Z2=A23
c let's get the LU decomposition of A11
call DGETRF( N2, N2, A22C, N2, IPIVZ2, INFO)
c let's solve the system A22*Z2=A23 to find Z2.
call DGETRS(TRANSN,N2,N3,A22C,N2,IPIVZ2,Z2,N2,INFO)
C=====
c we want to fins S, That's why we have to make some products:
c A31Z1=A31*Z1
call prod(N3,N1,N3,A31,Z1,A31Z1)
c A32Z2=A32*Z2
call prod(N3,N2,N3,A32,Z2,A32Z2)
c and now we can calculate S as:
do i=1,N3
do j=1,N3
S(i,j)=A33(i,j)-A31Z1(i,j)-A32Z2(i,j)
enddo
enddo
C===============================
C we need also to find the products: Z3=A31*(A11-1) and Z4=A32*(A22-1)
c we know that Z3*A11=A31===>(A11)'*(Z3)'=(A31)'
c the first thing we have to do is transposing A31:
call transp(N3,N1,A31,A31T)
c and we solve the system making the LU decomposition before.
call DGETRS(TRANST,N1,N3,A11C,N1,IPIVZ1,A31T,N1,INFO)
c now, A31T is Z3 trasposed. Let's find Z3:
call transp(N1,N3,A31T,Z3)
c====================
c let's do the same with Z4
call transp(N3,N2,A32,A32T)
c and we solve the system making the LU decomposition before.
call DGETRS(TRANST,N2,N3,A22C,N2,IPIVZ2,A32T,N2,INFO)
c now, A32T is Z4 trasposed. Let's find Z4:
call transp(N2,N3,A32T,Z4)
C we have all the elements of the LU decomposition.
end

C===========================================

subroutine solve(A11,A13,A22,A23,A31,A32,A33,
C N1,N2,N3,S,Z3,Z4,B1,B2,B3,A11C,A22C,IPIVZ1,IPIVZ2)
integer INFO,i,j,k
CHARACTER(1) TRANSN
integer IPIVZ3(N3),IPIVZ1(N1),IPIVZ2(N2)
double precision A11,A13,A22,A23,A31,A32,A33,A,A11C,A22C
C A31Z1,A32Z2,S,A31T,A32T,Z3,Z4,Z4B2,Z3B1,B1,B2,B3,A13X3,
C A23X3
dimension A11(N1,N1),A13(N1,N3),A22(N2,N2),A23(N2,N3),
C A31(N3,N1),A32(N3,N2),A33(N3,N3),A(N1+N2+N3,N1+N2+N3),
C A31Z1(N3,N3),A32Z2(N3,N3),S(N3,N3),Z3(N3,N1),Z4(N3,N2),
C Z4B2(N3),Z3B1(N3),B1(N1),B2(N2),B3(N3),A13X3(N1),A23X3(N2),
c A11C(N1,N1),A22C(N2,N2)
TRANSN='N'
N=N1+N2+N3
c let's find the solutions:
c we know that S*X3=B3-Z3*B1-Z4*B2
c we have to do the LU decomposition of S
c then we have to solve the system. Be compute B3=B3-Z3*B1-Z4*B2 and
c after calling DGETRS, B3 will be the solution X3.
call prod(N3,N1,1,Z3,B1,Z3B1)
call prod(N3,N2,1,Z4,B2,Z4B2)
do i=1,N3
B3(i)=B3(i)-Z3B1(i)-Z4B2(i)
enddo
call DGETRF(N3,N3,S,N3,IPIVZ3,INFO)
call DGETRS(TRANSN,N3,1,S,N3,IPIVZ3,B3,N3,INFO)
c now, we have the solutions of X3 in B3
c to calculate x1, we can compute: A11*X1=B1-A13*X3
call prod(N1,N3,1,A13,B3,A13X3)
do i=1,N1
B1(i)=B1(i)-A13X3(i)
enddo
call DGETRS(TRANSN,N1,1,A11C,N1,IPIVZ1,B1,N1,INFO)
c now, the solution X1 is in B1
c to calculate x2, we can copute: A22*X2=B2-A23*X3
call prod(N2,N3,1,A23,B3,A23X3)
do i=1,N2
B2(i)=B2(i)-A23X3(i)
enddo
call DGETRS(TRANSN,N2,1,A22C,N2,IPIVZ2,B2,N2,INFO)
c now, the solution X2 is in B2
end
C======================
C======================
SUBROUTINE prod(M1,M2,M3,A,B,C)
double precision A,B,C
integer M1,M2,M3
dimension A(M1,M2),B(M2,M3),C(M1,M3)
c let's comput C=A*B
do i=1,M1
do j=1,M3
C(i,j)=0
do k=1,M2
C(i,j)=C(i,j)+A(i,k)*B(k,j)
enddo
enddo
enddo
end
C======================
C======================
SUBROUTINE transp(M1,M2,A,B)
double precision A,B
integer M1,M2
dimension A(M1,M2),B(M2,M1)
c let's transpose A and put the values in B
do i=1,M2
do j=1,M1
B(i,j)=A(j,i)
enddo
enddo
end