\documentstyle[12pt,psfig]{article}
\topmargin=-0.2in\oddsidemargin=0in\textheight=8.8in\textwidth=6.5in  %format for exams
\newcommand{\bfl}{\begin{flushleft}}
\newcommand{\efl}{\end{flushleft}}
\newcommand{\ds}{\displaystyle}

\title{{ {\bf Software for the Adaptive Solution of Fredholm Integral Equations of the Second Kind (Regular and Singular Kernels)}}}
\author{{\sf {\bf Beny Neta}}\\
{\sf Naval Postgraduate School}\\ 
{\sf Department of Mathematics}\\
{\sf Code MA/Nd}\\ 
{\sf Monterey, CA 93943}}
\date{\today}

\begin{document}

\maketitle

\newpage

\begin{center}
\noindent {\bf Abstract}
\end{center}

Software for the adaptive solution of Fredholm integral equations of the second kind is given. The software allows the solution of equations having regular and singular kernels. The methods were developed by Neta and Nelson (1987) for regular kernels and Neta (1984) for singular ones. We discuss two types of singular kernels, i.e.  $\ds\log\,\vert\,x-y\,\vert$ and $\ds \vert\,x-y\,\vert^{-\alpha},\,\alpha\,<\,1.$  The choice of mesh points is made automatically so as to equidistribute both the change in the discrete solution and its gradient.  

\newpage
\section{\bf Regular Kernels}
\vskip 24pt
This section provides a FORTRAN program to solve Fredholm integral equations of the second kind having regular kernels. We give the main program, all subroutines required, the makefile and an example of input and output file. The solution is plotted via Matlab. We also plot the location of grid points in each stage of the process.

\bigskip
\begin{verbatim}
c
c     fiesk.f
c
c   Fredholm integral equation of second kind
c           regular kernel
c
c      need to supply kernel (rkf.f),
c                     first derivative of kernel (rkfx.f)
c                     second derivative of kernel (rkfxx.f)
c                     right hand side (gg.f)
c                     first derivative of rhs (gx.f)
c                     second derivative of rhs (gxx.f)
c                     exact solution (only for comparison) (exsol.f)
c
       real x(100),f(100),g(100),rk(100,100),d(100),fo(100)
       real rkd(100,100),rkdf(100),rkdfo(100),fxi(100),fxxi(100)
       real rkx(100),gsv(100),exct(100)
c
c     maxn    maximum number of nodes allowed (see dimension)
c     gamma   used in finding intervals where integral over
c             (x(j),x(j+1)) of absolute f prime is larger
c             than gamma * max abs f
c     delta   used in finding intervals where integral over
c             (x(j),x(j+1)) of absolute f double prime is larger
c             than delta * max abs f prime
c     c       the ratio between the largest and smallest interval
c
c
       read (5,*) maxn,gamma,delta,c
c2      format(i3,2f7.6,f3.1)
c
c     indprt = 1  print all
c            = 0  print some
c     niter   max number of iterations used by Gauss Seidel
c     tol     tolerance for convergence
c
        read(5,*) indprt,niter,tol
c1       format(i1,i3,f15.10)
c
c     a,b   endpoints of integration interval
c
        read(5,*) a,b
c10      format(2f10.5)
c
c     n    number of nodes to start with (default=9)
c     m    number of additional nodes (default=91)
c     maxal maximum nodes allowed (n+m)
c
        read(5,*,end=30) n,m
c20      format(2i2)
        go to 40
30      n=9
        m=maxn-n
40      continue
        if(n.gt.maxn) n=maxn
        if((n+m).gt.maxn) m=maxn-n
        maxal=n+m
        write(6,50) n,a,b,m
50      format('1'//10x,'Approximate a Fredholm integral equation',
     &   ' of the second kind using ',i2,' equally spaced nodes'//
     &    15x,'the interval of integration is (',f10.5,',',f10.5,
     &    ')'//15x,'the program uses adaptive trapezoidal rule'//
     &  15x,'dont use more than ',i2,' extra nodes'/)
c
c  *** THIS FORMAT SHOULD BE CHANGED WHEN KERNEL, RIGHT HAND SIDE,
c      SOLUTION ARE CHANGED
c
        write(6,52)
52      format(15x,'The kernel is          k(x,y)=(x-.5)*(y-.5)'//
     &         15x,'the right hand side is g(x)=exp(-10*(x-.5)**2)'//
     &         15x,'the exact solution is  f(x)=exp(-10*(x-.5)**2)'/)
        write(6,54) maxn,gamma,delta,c
54      format(10x,'maximum number of nodes allowed (dimension) is',
     &   i4/10x,'the constant gamma used in finding intervals where',
     &   ' integral over abs f prime is large ',f14.7/10x,
     &   'the constant delta used in finding intervals where',
     &   ' integral over abs f double prime is large ',f14.7/10x,
     &    'the ratio between the largest and smallest interval ',f14.7/)
        do 55 i=1,n
        f(i)=0.
        fo(i)=0.
55      continue
c
c   h  length of initial subintervals
c
        h=(b-a)/(n-1.)
c
c   nodes are equally spaced
c
        x(1)=a
        do 60 i=2,n
60      x(i)=x(i-1)+h
        if (indprt.eq.1) print 62,(x(i),i=1,n)
62      format(2x,'x ',5e14.7)
        write(30,*) n
        write(30,*) (x(i),i=1,n)
c
c     compute right hand side
c
        call gfcn(x,n,g,0)
        if (indprt.eq.1) print 64,(g(i),i=1,n)
64      format(2x,'g ',5e14.7)
c
c     compute exact solution
c
        call exsol(x,n,exct)
        if (indprt.eq.1) print 640,(exct(i),i=1,n)
640      format(2x,'exact ',5e14.7)
c
c    initialize rkdfo = fo - g (the integral term)
c
        do 65 i=1,n
65      rkdfo(i)=fo(i)-g(i)
c
c    compute k(x,y)
c
        call kfcn(x,n,rk,maxn,0)
        if (indprt.eq.1) print 66,((rk(i,j),i=1,n),j=1,n)
66      format(2x,'rk ',5e14.7)
c
c   d gives the weights
c
        call dfcn(x,n,a,b,d,0)
        if (indprt.eq.1) print 68,(d(i),i=1,n)
68      format(2x,'d ',5e14.7)
c
c    compute the entries of the matrix
c
        call rkdfcn(rk,d,n,maxn,rkd)
        if (indprt.eq.1) print 70,((rkd(i,j),i=1,n),j=1,n)
70      format(2x,'rkd ',5e14.7)
c
c   Gauss Seidel to solve the system
c
        call gsitr(rkd,g,f,n,maxn,niter,gsv,tol)
        write(6,230)
        do 80 i=1,n
        xx=x(i)
        if (i.lt.n) h=x(i+1)-xx
        diff=f(i)-exct(i)
        write(6,220) xx,exct(i),f(i),diff,h
80      continue
c
c   evaluate rkdf (the integral term)
c
        do 82 i=1,n
82      rkdf(i)=f(i)-g(i)
        if(m.lt.1) go to 200
c
c   loop to add up to m nodes when and where necessary
c
        do 100 kk=1,m
c
c   find largest f in abs value. denote it fl
c
        fl=abs(f(1))
        do 69 l=2,n
        if(fl.lt.abs(f(l))) fl=abs(f(l))
69      continue
        fl=delta*fl
        nm1=n-1
c
c   evaluate fxi(j) j=2,...n-1
c   which is the integral of absolute value of f prime
c   over (x(j),x(j+1)).
c
        fxl=0.
        do 84 j=1,nm1
        call kxfcn(x,j,n,rkx)
        x1=0.
        x2=0.
        do 86 l=1,nm1
        x1=x1+(rkx(l)*f(l)+rkx(l+1)*f(l+1))*(x(l+1)-x(l))
86      continue
        x1=x1*.5+gx(x(j))
c
c   compute absolute f prime
c
        if(fxl.lt.abs(x1)) fxl=abs(x1)
        jp1=j+1
        call kxfcn(x,jp1,n,rkx)
        do 88 l=1,nm1
        x2=x2+(rkx(l)*f(l)+rkx(l+1)*f(l+1))*(x(l+1)-x(l))
88      continue
        x2=x2*.5+gx(x(jp1))
        fxi(j)=abs(x1+x2)*.5*(x(jp1)-x(j))
84      continue
        if (indprt.eq.1) print 300,(fxi(j),j=1,nm1)
300     format(2x,'fxi',5e14.7)
        fxl=fxl*gamma
c
c   evaluate fxxi(j) j=2,...n-1
c   which is the integral of absolute value of f double prime
c   over (x(j),x(j+1)).
c
        do 74 j=1,nm1
        call kxxfcn(x,j,n,rkx)
        x1=0.
        x2=0.
        do 76 l=1,nm1
        x1=x1+(rkx(l)*f(l)+rkx(l+1)*f(l+1))*(x(l+1)-x(l))
76      continue
        x1=x1*.5+gxx(x(j))
        jp1=j+1
        call kxxfcn(x,jp1,n,rkx)
        do 78 l=1,nm1
        x2=x2+(rkx(l)*f(l)+rkx(l+1)*f(l+1))*(x(l+1)-x(l))
78      continue
        x2=x2*.5+gxx(x(jp1))
        fxxi(j)=abs(x1+x2)*.5*(x(jp1)-x(j))
74      continue
        if (indprt.eq.1) print 302,(fxxi(j),j=1,nm1)
302     format(2x,'fxxi',5e14.7)
c
c   compute the max norm of difference between 2 consecutive
c   approximate solutions.
c  
c   rinf   norm
c   indx   index where it occures
c
        call infnrm(f,fo,rkdf,rkdfo,n,rinf,indx,rkinf,kindx)
106     continue
c
c   find the smallest interval length
c
        hmin=h
        do 94 j=1,nm1
        if((x(j+1)-x(j)).lt.hmin) hmin=x(j+1)-x(j)
94      continue
        kl=0
        jj=0
        do 102 j=1,nm1
c
c   any interval larger than c*hmin should be cut
c
        if((x(j+jj+1)-x(j+jj)).gt.c*hmin) go to 96
        if(fxi(j).lt.fl.and.fxxi(j).lt.fxl) go to 102
96      kl=1
        indx=j+jj+1
        if(n.ge.maxal) go to 103
c
c   shift routine moves entries of f,g,x,k to make room for the 
c   additional node
c
         call shift(x,f,fo,d,g,rk,n,maxn,indx)
        jj=jj+1
        if (indprt.eq.1) print 62,(x(i),i=1,n)
c
c  compute right hand side for the additional node
c
        call gfcn(x,n,g,indx)
        if (indprt.eq.1) print 64,(g(i),i=1,n)
c
c  compute additional rows and columns of k(x,y) for the additional node
c
        call kfcn(x,n,rk,maxn,indx)
        if (indprt.eq.1) print 66,((rk(i,ij),i=1,n),ij=1,n)
c
c  compute weights for additional node
c
        call dfcn(x,n,a,b,d,indx)
        if (indprt.eq.1) print 68,(d(i),i=1,n)
c
c  compute entries of matrix
c
        call rkdfcn(rk,d,n,maxn,rkd)
        if (indprt.eq.1) print 70,((rkd(i,ij),i=1,n),ij=1,n)
102     continue
103     if(kl.eq.0) go to 200
c
c   update the exact solution
c
       call exsol(x,n,exct)
       if (indprt.eq.1)  print 640,(exct(i),i=1,n)
c
c   solve by Gauss Seidel the new system using the previous solution
c   as initial guess. At the new node, the guess is the average of
c   solutions at neighboring points
c
        call gsitr(rkd,g,f,n,maxn,niter,gsv,tol)
c
c   print exact solution, approximate solution and difference
c
        write(6,230)
230     format(8x,'node',7x,'exact solution  approx.solution   ',
     &   'difference'14x,'h'/)
        do 240 i=1,n
        xx=x(i)
        if(i.lt.n) h=x(i+1)-xx
        diff=f(i)-exct(i)
        write(6,220) xx,exct(i),f(i),diff,h
220     format(3x,e14.7,2x,e14.7,2x,e14.7,2x,e14.7,5x,e14.7)
240     continue
        write(30,*) n
        write(30,*) (x(i),i=1,n)
c
c   evaluate rkd (the integral term)
c
        do 92 i=1,n
92      rkdf(i)=f(i)-g(i)
        if(n.ge.maxal) go to 104
100     continue
c
c   process did not converge using n+m nodes
c
104     continue
c
c   compute the max norm of difference between 2 consecutive
c   approximate solutions and it's location
c  
        call infnrm(f,fo,rkdf,rkdfo,n,rinf,indx,rkinf,kindx)
        write(6,110) n,rinf,indx
110     format(10x,'no convergence using ',i3,' nodes'/10x,
     &   'the max error is ',e14.7,' at component ',i3,
     &   ' of vector f')
200     continue
        stop
        end
\end{verbatim}

\bigskip
\begin{verbatim}
       subroutine dfcn(x,n,a,b,d,ind)
       real x(1),d(1)
c
c   subroutine to compute weights
c
       if(ind.gt.0) go to 150
       d(1)=(x(2)-x(1))*(b-a)*.5
       d(n)=(x(n)-x(n-1))*(b-a)*.5
       nm1=n-1
       do 100 i=2,nm1
100    d(i)=(x(i+1)-x(i-1))*(b-a)*.5
       return
150    continue
       d(ind-1)=(x(ind  )-x(ind-2))*(b-a)*.5
       d(ind  )=(x(ind+1)-x(ind-1))*(b-a)*.5
       d(ind+1)=(x(ind+2)-x(ind  ))*(b-a)*.5
       if(ind.eq.2) d(ind-1)=(x(ind  )-x(ind-1))*(b-a)*.5
       if(ind.eq.nm1) d(ind+1)=(x(ind+1)-x(ind))*(b-a)*.5
       return
       end
\end{verbatim}

\bigskip
\begin{verbatim}
       subroutine exsol(x,n,exct)
       real x(1),exct(1)
       ad1=exp(1./12.)
       ad=3.*(ad1-1.)
       a=1./ad
       b=(ad1-2.)/ad
       cd1=exp(3./8.)
       cd2=exp(.25)
       c=13./(3.*(cd1-cd2))
       d=(2.*cd1-15.*cd2)/(3.*(cd1-cd2))
       do 100 i=1,n
       xx=x(i)
       if(xx.ge.0. .and. xx.le. .25) exct(i)=1./3.
       if(xx.gt..25 .and. xx.lt. (1./3.)) exct(i)=a*exp(xx-.25)+b
       if(xx.ge.(1./3.) .and. xx.le. .5) exct(i)=2./3.
       if(xx.gt..5 .and. xx.lt. (5./8.)) exct(i)=c*exp(xx-.25)+d
       if(xx.ge.(5./8.) .and. xx.le. 1.) exct(i)=5.
100    continue
       return
       end
\end{verbatim}

\bigskip
\begin{verbatim}
       subroutine gfcn(x,n,g,ind)
       real x(1),g(1)
c
c   subroutine to compute right hand side
c
       if(ind.gt.0) go to 150
       do 100 i=1,n
100    g(i)=gg(x(i))
       return
150    continue
       g(ind-1)=gg(x(ind-1))
       g(ind  )=gg(x(ind  ))
       g(ind+1)=gg(x(ind+1))
       return
       end
\end{verbatim}

\bigskip
\begin{verbatim}
       function gg(xx)
       ad1=exp(1./12.)
       ad=3.*(ad1-1.)
       a=1./ad
       b=(ad1-2.)/ad
       cd1=exp(3./8.)
       cd2=exp(.25)
       cd=3.*(cd1-cd2)
       c=13./cd
       d=(2.*cd1-15.*cd2)/cd
       cons=5461./3456.-(185.*ad1-216.)/(864.*(ad1-1.))-
     &        (606.*cd1-697.*cd2)/(384.*(cd1-cd2))
       if(xx.ge.0. .and. xx.le. .25) exc=1./3.
       if(xx.gt..25 .and. xx.lt. (1./3.)) exc=a*exp(xx-.25)+b
       if(xx.ge.(1./3.) .and. xx.le. .5) exc=2./3.
       if(xx.gt..5 .and. xx.lt. (5./8.)) exc=c*exp(xx-.25)+d
       if(xx.ge.(5./8.) .and. xx.le. 1.) exc=5.
       gg=exc-cons*xx**2
       return
       end
\end{verbatim}

\bigskip
\begin{verbatim}
       subroutine gsitr(a,b,x,n,ndim,niter,g,tol)
       real x(1),b(1),g(1),a(ndim,ndim)
c
c  Gauss Seidel iterative method to solve the system
c
       do 1 i=1,n
1      g(i)=b(i)
       do 10 i=1,n
       save=a(i,i)
       b(i)=b(i)/save
       do 5 j=1,n
       a(i,j)=a(i,j)/save
5      continue
10     continue
       do 40 iter=1,niter
       xmax=0.
       do 30 i=1,n
       save=x(i)
       x(i)=b(i)
       do 20 j=1,n
       if(j.eq.i) go to 20
       x(i)=x(i)-a(i,j)*x(j)
20     continue
       if(abs(x(i)-save).gt.xmax) xmax=abs(x(i)-save)
30     continue
       if(xmax.le.tol) go to 300
40     continue
       write(6,200) tol,niter
200    format('0','did not meet tolerance of ',e14.7,' in ',
     &   i4,' iterations'/1x,'last value of x was returned to',
     &   ' caller')
300    continue
       print 400,iter,n
400    format(30x,'number of iterations used by Gauss Seidel is ',i4,
     &   10x,'n=',i4)
       do 500 i=1,n
500    b(i)=g(i)
       return
       end
\end{verbatim}

\bigskip
\begin{verbatim}
       function gx(xx)
       ad1=exp(1./12.)
       ad=3.*(ad1-1.)
       a=1./ad
       b=(ad1-2.)/ad
       cd1=exp(3./8.)
       cd2=exp(.25)
       cd=3.*(cd1-cd2)
       c=13./cd
       d=(2.*cd1-15.*cd2)/cd
       cons=5461./3456.-(185.*ad1-216.)/(864.*(ad1-1.))-
     &        (606.*cd1-697.*cd2)/(384.*(cd1-cd2))
       if(xx.ge.0. .and. xx.le. .25) exc=0.
       if(xx.gt..25 .and. xx.lt. (1./3.)) exc=a*exp(xx-.25)+b
       if(xx.ge.(1./3.) .and. xx.le. .5) exc=0.
       if(xx.gt..5 .and. xx.lt. (5./8.)) exc=c*exp(xx-.25)+d
       if(xx.ge.(5./8.) .and. xx.le. 1.) exc=0.
       gx=exc-2.*cons*xx
       return
       end
\end{verbatim}

\bigskip
\begin{verbatim}
       function gxx(xx)
       ad1=exp(1./12.)
       ad=3.*(ad1-1.)
       a=1./ad
       b=(ad1-2.)/ad
       cd1=exp(3./8.)
       cd2=exp(.25)
       cd=3.*(cd1-cd2)
       c=13./cd
       d=(2.*cd1-15.*cd2)/cd
       cons=5461./3456.-(185.*ad1-216.)/(864.*(ad1-1.))-
     &        (606.*cd1-697.*cd2)/(384.*(cd1-cd2))
       if(xx.ge.0. .and. xx.le. .25) exc=0.
       if(xx.gt..25 .and. xx.lt. (1./3.)) exc=a*exp(xx-.25)
       if(xx.ge.(1./3.) .and. xx.le. .5) exc=0.
       if(xx.gt..5 .and. xx.lt. (5./8.)) exc=c*exp(xx-.25)
       if(xx.ge.(5./8.) .and. xx.le. 1.) exc=0.
        gxx=exc-2.*cons
       return
       end
\end{verbatim}

\bigskip
\begin{verbatim}
       subroutine infnrm(f,fo,rkdf,rkdfo,n,rinf,indx,rkinf,kindx)
       real f(1),fo(1),rkdf(1),rkdfo(1)
c
c  subroutine to compute the max difference between two consecutive
c  approximations and its location
c
       rinf=0.
       indx=1
       rkinf=0.
       kindx=1
       do 100 i=1,n
       qk=abs(rkdf(i)-rkdfo(i))
       q=abs(f(i)-fo(i))
       rkdfo(i)=rkdf(i)
       fo(i)=f(i)
       if(rkinf.gt.qk) go to 50
       rkinf=qk
       kindx=i
50     continue
       if(rinf.gt.q) go to 100
       rinf=q
       indx=i
100    continue
       return
       end
\end{verbatim}

\bigskip
\begin{verbatim}
       subroutine kfcn(x,n,rk,ndim,ind)
       real x(1),rk(ndim,ndim)
c
c  subroutine to compute k(x,y)
c
       if(ind.gt.0) go to 150
       do 100 i=1,n
       xx=x(i)
       do 100 j=1,n
       yy=x(j)
100    rk(i,j)=rkf(xx,yy)
       return
150    continue
       indm1=ind-1
       indp1=ind+1
       do 200 i=1,n
       xx=x(i)
       do 200 j=indm1,indp1
       yy=x(j)
       rk(j,i)=rkf(yy,xx)
       if(i.lt.indm1.or.i.gt.indp1) rk(i,j)=rkf(xx,yy)
200    continue
       return
       end
\end{verbatim}

\bigskip
\begin{verbatim}
       subroutine kxfcn(x,j,n,rkx)
       real x(1),rkx(1)
c
c  subroutine to compute first derivative of k with respect to x
c
       xj=x(j)
       do 100 i=1,n
       xi=x(i)
100    rkx(i)=rkfx(xj,xi)
       return
       end
\end{verbatim}

\bigskip
\begin{verbatim}
       subroutine kxxfcn(x,j,n,rkx)
       real x(1),rkx(1)
c
c  subroutine to compute second derivative of k with respect to x
c
       xj=x(j)
       do 100 i=1,n
       xi=x(i)
100    rkx(i)=rkfxx(xj,xi)
       return
       end
\end{verbatim}

\bigskip
\begin{verbatim}
       subroutine rkdfcn(rk,d,n,ndim,rkd)
       real d(1),rk(ndim,ndim),rkd(ndim,ndim)
c
c  subroutine to compute entries of the matrix
c
       do 100 i=1,n
       dd=d(i)
       do 100 j=1,n
       rkd(j,i)=-rk(j,i)*dd
       if(i.eq.j)  rkd(i,i)=rkd(i,i)+1.
100    continue
       return
       end
\end{verbatim}

\bigskip
\begin{verbatim}
       function rkf(xx,yy)
       rkf=xx**2*yy
       return
       end
\end{verbatim}

\bigskip
\begin{verbatim}
       function rkfx(xx,yy)
       rkfx=2.*xx*yy
       return
       end
\end{verbatim}

\bigskip
\begin{verbatim}
       function rkfxx(xx,yy)
       rkfxx=2.*yy
       return
       end
\end{verbatim}

\bigskip
\begin{verbatim}
       subroutine shift(x,f,fo,d,g,rk,n,ndim,indx)
       real x(1),f(1),fo(1),d(1),g(1),rk(ndim,ndim)
c
c  subroutine to shift entries of x,f,d,g,k(x,y) to make room
c   for additional node
c
       np1=n+1
       do 100 i=indx,n
       nn=n+indx-i
       x(nn+1)=x(nn)
       f(nn+1)=f(nn)
       fo(nn+1)=fo(nn)
       d(nn+1)=d(nn)
       g(nn+1)=g(nn)
       do 100 j=1,n
       rk(j,nn+1)=rk(j,nn)
100    continue
       do 200 i=indx,n
       nn=n+indx-i
       do 200 j=1,np1
       rk(nn+1,j)=rk(nn,j)
200    continue
       x(indx)=(x(indx-1)+x(indx+1))*.5
       f(indx)=(f(indx-1)+f(indx+1))*.5
       fo(indx)=f(indx)
       n=np1
       return
       end
\end{verbatim}

\bigskip
We now give an example of input file.
The first line gives:
\begin{itemize}
\item maxn=maximum number of nodes allowed (see dimension)
\item  gamma   used in finding intervals where integral over
$(x(j),x(j+1))$ of $\vert f' \vert$ is larger
than gamma * max $\vert f \vert$
\item delta   used in finding intervals where integral over
$(x(j),x(j+1))$ of $\vert f'' \vert$ is larger
than delta * max $\vert f' \vert$
\item c       the ratio between the largest and smallest interval
\end{itemize}

\medskip
The second line gives:
\begin{itemize}
\item indprt = 1  print all
indprt = 0  print some

\item niter   max number of iterations used by Gauss Seidel
\item tol     tolerance for convergence

\end{itemize}

\medskip
The third line gives:
a,b   endpoints of integration interval

\medskip
The fourth line gives:
\begin{itemize}
\item  n    number of nodes to start with (default=9)
\item m    number of additional nodes (default=91)
\end{itemize}

\medskip
\begin{verbatim}
100 .1 .1 10.
0 100 .00000001
0. 1.
5 95
\end{verbatim}

\bigskip
Here is the makefile used
\medskip
\begin{verbatim}
OBJ= fiesk.o  dfcn.o gfcn.o gg.o gsitr.o gx.o gxx.o infnrm.o \
                     kfcn.o kxfcn.o kxxfcn.o rkdfcn.o \
                     rkf.o rkfx.o rkfxx.o shift.o exsol.o
run_fiesk : $(OBJ)
	f77 -C -O -r8 -o $@ $(OBJ)
\end{verbatim}

\bigskip
We now give the Matlab program used to plot the solution. The program reads one of the output files (fort.30) and extracts the necessary data. The first figure is the location of the grid points at each stage (assuming 6 stages were used), The second figure is the approximate solution at the last stage and the last figure is the exact solution.
\medskip

\begin{verbatim}
clg
grdx=[0.,.2,.4,.6,.8,1.];
grdy=[-1.,0.,1.,2.];
       ad1=exp(1./12.);
       ad=3.*(ad1-1.);
       a=1./ad;
       b=(ad1-2.)/ad;
       cd1=exp(3./8.);
       cd2=exp(.25);
       c=13./(3.*(cd1-cd2));
       d=(2.*cd1-15.*cd2)/(3.*(cd1-cd2));
fid=fopen('fort.30','r');
ab=fscanf(fid,'%d',1);
n=ab(1)
[aa,count]=fscanf(fid,'%e%e%e%e%e',n);
ab(1:n)=ones(1,n);
plot(aa,ab,'+')
hold on
x1=aa;
n1=n;
clear aa;
clear ab;
ab=fscanf(fid,'%d',1);
n=ab(1)
[aa,count]=fscanf(fid,'%e%e%e%e%e',n);
ab(1:n)=1.1*ones(1,n);
plot(aa,ab,'+')
x2=aa;
n2=n;
clear aa;
clear ab;
ab=fscanf(fid,'%d',1);
n=ab(1)
[aa,count]=fscanf(fid,'%e%e%e%e%e',n);
ab(1:n)=1.2*ones(1,n);
plot(aa,ab,'+')
x3=aa;
n3=n;
clear aa;
clear ab;
ab=fscanf(fid,'%d',1);
n=ab(1)
[aa,count]=fscanf(fid,'%e%e%e%e%e',n);
ab(1:n)=1.3*ones(1,n);
plot(aa,ab,'+')
x4=aa;
n4=n;
clear aa;
clear ab;
ab=fscanf(fid,'%d',1);
n=ab(1)
[aa,count]=fscanf(fid,'%e%e%e%e%e',n);
ab(1:n)=1.4*ones(1,n);
plot(aa,ab,'+')
x5=aa;
n5=n;
clear aa;
clear ab;
ab=fscanf(fid,'%d',1);
n=ab(1)
[aa,count]=fscanf(fid,'%e%e%e%e%e',n);
ab(1:n)=1.5*ones(1,n);
plot(aa,ab,'+')
hold off
axis([0,1.,.1,3])
set(gca,'XTick',[grdx]);
%set(gca,'YTick',[grdy]);
print fig1 -dps
clear exct;
for i=1:n1
       xx=x1(i);
       if( xx <= .25) exct(i)=1./3.;
       elseif(xx < (1./3.)) exct(i)=a*exp(xx-.25)+b;
       elseif( xx <= .5) exct(i)=2./3.;
       elseif( xx < (5./8.)) exct(i)=c*exp(xx-.25)+d;
       else exct(i)=5.;
end
end
plot(x1,exct)
print fig2 -dps
clear exct;
clear x1;
for i=1:10*n5
       x1(i)=(i-1)/(10*n5-1);
       xx=x1(i);
       if( xx <= .25) exct(i)=1./3.;
       elseif(xx < (1./3.)) exct(i)=a*exp(xx-.25)+b;
       elseif( xx <= .5) exct(i)=2./3.;
       elseif( xx < (5./8.)) exct(i)=c*exp(xx-.25)+d;
       else exct(i)=5.;
end
end
plot(x1,exct)
print fig3 -dps
\end{verbatim}

\section{Example of Plots}

In this section we bring the three figures plotted by Matlab for the case described in this software.

\begin{figure}[htb]
\centerline{ \psfig{figure=fig1.ps,width=3.5in} }
\caption{The Location of Grid Points.}
\label{fig1}
\end{figure}

\newpage
\begin{figure}[htb]
\centerline{ \psfig{figure=fig2.ps,width=3.5in} }
\caption{The Approximate Solution at the Last Stage.}
\label{fig2}
\end{figure}

%\newpage
\begin{figure}[bth]
\centerline{ \psfig{figure=fig3.ps,width=3.5in} }
\caption{The Exact Solution.}
\label{fig3}
\end{figure}

\newpage
%\vskip 24pt

\newpage
\section{\bf Singular Kernels}
\vskip 24pt
This section provides a FORTRAN program to solve Fredholm integral equations of the second kind having singular kernels. Two types of singularity are considered, i.e.  $\ds\log\,\vert\,x-y\,\vert$ and $\ds \vert\,x-y\,\vert^{-\alpha},\,\alpha\,<\,1.$ We give the main program, all subroutines required, the makefile and an example of input and output file. 

\bigskip
\begin{verbatim}
\end{verbatim}

\bigskip
\begin{verbatim}
\end{verbatim}

\bigskip
\begin{verbatim}
\end{verbatim}

\bigskip
\begin{verbatim}
\end{verbatim}

\bigskip
\begin{verbatim}
\end{verbatim}

\bigskip
\begin{verbatim}
\end{verbatim}

\bigskip
\begin{verbatim}
\end{verbatim}

\bigskip
\begin{verbatim}
\end{verbatim}

\bigskip
\begin{verbatim}
\end{verbatim}

\bigskip
\begin{verbatim}
\end{verbatim}

\bigskip
\begin{verbatim}
\end{verbatim}


\newpage
%\vskip 24pt
\noindent {\bf References}

\vskip 12pt
\noindent 

\newpage
\begin{center}\bf Distribution List \rm\end{center}
\medskip \hfill No. of copies

\bfl
Director \hfill 2 ~~~~~~~~~ \\ 
Defense Technical Information Center\\
8725 John J. Kingman Rd., Ste 0944\\
Ft. Belvoir, VA \, 22060-6218
\efl

\medskip
\begin{flushleft}
Dean of Research  \hfill 1 ~~~~~~~~~ \\
Code 09\\
Naval Postgraduate School\\
Monterey, CA\, 93943
\end{flushleft}


\medskip
\begin{flushleft}
Library \hfill 2 ~~~~~~~~~ \\
Code 013\\
Naval Postgraduate School\\
Monterey, CA \, 93943
\end{flushleft}

\medskip
\begin{flushleft}
Professor Beny Neta \hfill 5 ~~~~~~~~~ \\
Code MA/Nd\\
Naval  Postgraduate School\\
Monterey, CA \, 93943
\end{flushleft}


\end{document}
