radial_basis_function.f90

4总常用的插值程序的Fortran90源代码

program rbf
!
!*******************************************************************************
!
!! RBF implements the radial basis function method.
!
!
!  mq1dlsc - a multi-quadric radial basis function method
!            with least squares method.
!
!  version "c" - the alfa-coefficients and c_j are to be found
!
!           test 1: 1d function interpolation
!           test 2: 1d elliptic equation
!
!
  integer, parameter :: nmax = 501

  integer, parameter :: lwa = ( nmax * ( 3 * nmax + 13 ) / 2 + 1 )
  integer, parameter :: nxtm = 1001
!
  real a(nmax,nmax)
  real aaa(nmax,nmax)
  real adfdx
  real anal
  real asa(nmax,nmax)
  real b(nmax)
  real c(nmax)
  real ck
  real cmax
  real cmin
  character ( len = 30 ) comnd
  real cond
  real d2fdx2
  real d2fundx
  real dfdx
  real dfundx
  real dsdx
  real dx
  real dxfun
  real dxt
  real err
  real errdf
  real f(nmax)
  real ffun
  real fun
  real fvec(nmax)
  real g
  real gx
  real g2x
  real hxt
  integer i
  integer ifun
  integer in
  integer info
  integer io
  integer io2
  integer ios
  integer ipvt(nmax)
  integer j
  integer k
  integer m
  integer mode
  integer n
  integer nn
  integer nnxt
  integer npmx
  integer nqq
  integer nxt
  character ( len = 10 ) par
  real radc
  real rcur
  real sol
  real tol
  real uxx
  real w(nmax)
  real wa(lwa)
  real xk(nmax)
  real x1
  real x2
  real xdx
  real xfun
  real xt(nxtm)
  real xt1
  real xt2
  real xxf
  real y(nmax)
  real z
!
  external hlsrec
!
  common /i_common/ ifun, m, mode, nn, nqq, nxt
  common /r_common/ aaa, asa, b, hxt, xk, xt
!
!  declaration operator-functions
!
  fun(z) = ffun(ifun,z)
  dfdx(z) = dfundx(ifun,z)
  d2fdx2(z) = d2fundx(ifun,z)
!
  g(j,z) = sqrt ( ( z - xk(j) )**2 + c(j)**2 )

  gx(j,z) = ( z - xk(j) ) / sqrt ( ( z - xk(j) )**2 + c(j)**2 )

  g2x(j,z) = 1.0 / sqrt ( ( z - xk(j) )**2 + c(j)**2 ) &
     - ( z - xk(j) )**2 / sqrt ( ( ( z - xk(j) )**2 + c(j)**2 )**3 )
!
!  Read the solver parameters.
!
  ifun = 1
  in = 5
  io = 9
  io2 = 10
  open ( io, file='mqle.dat', status = 'unknown' )

  mode = 2

  write ( *, * ) ' '
  write ( *, * ) 'Choose the function:'
  write ( *, * ) '  1 = x**2'
  write ( *, * ) '  2 = x*(x-1)**2'
  write ( *, * ) '  3 = sin(2*pi*x)'
  write ( *, * ) '  4 = exp(-x)'
  write ( *, * ) '  5 = abs(x-0.5)'
  write ( *, * ) '  6 = atan(10*(x-0.5)): '
  read ( *, * ) ifun

  write ( *, * ) ' '
  write ( *, * ) 'Enter N:'
  read ( *, * ) n

  write ( *, * ) ' '
  write ( *, * ) 'Enter NNXT'
  read ( *, * ) nnxt

  nxt = abs ( nnxt )
  nxt = ( nxt - 1 ) / 2 * 2 + 1
  m = n

  write ( *, * ) ' '
  write ( *, * ) 'Enter K, CMIN, CMAX:'
  write ( *, * ) ' '
  write ( *, * ) '  K=1: C(J)=cmin*(j/n)^(3/4);'
  write ( *, * ) '  K=2: C(J)=cmin*((cmax/cmin)^((j-1)/(n-1)))'
  read ( *, * ) k, cmin, cmax
!
!  Generate CJ
!
  do j = 1, n
    if ( k == 1 ) then
      if ( cmax > 0.0 ) then
        c(j) = cmax * ( real ( j ) / real ( n ) )**0.75
      else
        c(j) = abs(cmax)
      end if
    else
      c(j) = cmin * ( cmax / cmin )**( float ( j - 1 ) / float ( n - 1 ) )	
    end if
  end do

  write ( *, * ) ' '
  write ( *, * ) 'Shape parameters C(1:N):'
  write ( *, * ) ' '
  write ( *, '(i3,g12.4)' ) ( i, c(i), i = 1, n )
!
!  Set integration points
!
  xt1 = 0.0
  xt2 = 1.0
  xt(1) = xt1
  dxt = ( xt2 - xt1 ) / real ( nxt - 1 )
  hxt = dxt
  do i = 2, nxt
    xt(i) = xt(i-1) + dxt
  end do

  write ( *, * ) ' '
  write ( *, * ) 'Integration points XT(1:NXT):'
  write ( *, * ) ' '
  write ( *, '(i3,g12.4)' ) ( i, xt(i), i = 1, nxt )
!
!  Set xj points - function knots
!
  x1 = 0
  x2 = 1
  xk(1) = x1
  dx = ( x2 - x1 ) / real ( n - 1 )
  do i = 2, n
    xk(i) = xk(i-1) + dx
  end do

  write ( *, * ) ' '
  write ( *, * ) 'Function points XK(1:N):'
  write ( *, * ) ' '
  write ( *, '(i3,g12.4)' ) ( i, xk(i), i = 1, n )
!
!  form a matrix
!
610  continue

  nn = n

  a(1:nn,1:nn) = 0.0

  do j = 1, n
    do i = 1, n
      a(i,j) = g(j,xk(i))
    end do
  end do
!
!  Form a right hand side.
!
  do i = 1, n
    b(i) = fun ( xk(i) )
  end do
!
!   stiff b.c. (function at the boundary)
!
  if ( nnxt < 0 ) then
!
!   set b.c.
!
    i = 1
    do j = 1, n
      a(i,j) = g(j,xk(i))
      asa(i,j) = a(i,j)
    end do
    b(i) = fun ( xk(i) )

    i = n
    do j = 1, n
      a(i,j) = g(j,xk(i))
      asa(i,j) = a(i,j)
    end do
    b(i) = fun ( xk(i) )

  end if
!
!  Decompose the system matrix.
!
  call decomp ( nmax, n, a, cond, ipvt, w )

  write ( *, * ) ' '
  write ( *, * ) '  Number of knots = ', n
  write ( *, * ) '  Number of points = ', nxt
  write ( *, * ) '  Condition number = ', cond
!
!  solve by hmin/hybrid
!
  y(1:n) = b(1:n)

  call solve ( nmax, n, a, y, ipvt )

  write ( *, * ) ' '
  write ( *, * ) 'Preliminary ALFA(1:N)'
  write ( *, * ) ' '
  write ( *, '(i3,e14.3)' ) ( i, y(i), i = 1, n )
  tol = 1.d-12
  info = 0
  y(1:n) = c(1:n)

  call hybrd1 ( hlsrec, nn, y, fvec, tol, info, wa, lwa )

  write(*,2000) tol, info
2000  format(' solved o.k.',' tol req =',1pe12.3,' info =',i5)

  if ( info == 0 ) then
    write ( *, * ) 'improper input parameters.'
  else if ( info == 1 ) then
    write ( *, * ) 'relative error between x and the solution is at most tol.'
  else if ( info == 2 ) then
    write ( *, * ) 'number of calls to fcn has reached or exceeded 200*(n+1).'
  else if ( info == 3 ) then
    write ( *, * ) 'tol is too small.'
    write ( *, * ) 'no further improvement in the solution x is possible.'
  else if ( info == 4 ) then
    write ( *, * ) 'iteration is not making good progress.'
  end if

  write ( *, * ) ' '
  write ( *, * ) 'Distance L(1:N):'
  write ( *, * ) ' '
  write ( *, '(i3,e14.3)' ) ( i, y(i), i = 1, n )
  read ( *, * )

  c(1:n) = y(1:n)
!
!  Solve again using new c
!
  do j = 1, n
    do i = 1, n
      a(i,j) = g(j,xk(i))
      asa(i,j) = a(i,j)
    end do
  end do

  do i = 1, n
    b(i) = fun ( xk(i) )
  end do

  call decomp ( nmax, n, a, cond, ipvt, w )

  write ( *, * ) ' '
  write ( *, * ) '  Number of knots = ', n
  write ( *, * ) '  Number of points = ', nxt
  write ( *, * ) '  Condition number = ', cond

  y(1:n) = b(1:n)

  call solve ( nmax, n, a, y, ipvt )

  write ( *, * ) ' '
  write ( *, * ) 'Final ALFA(1:N)'
  write ( *, * ) ' '
  write ( *, '(i3,e14.3)' ) ( i, y(i), i = 1, n )
  read ( *, * ) 
!
!  print solution at the nodes
!
!cc  write (*,3000)n

  rewind io
  write ( *, 3000 ) n, cond
3000  format('# solution ',i5,' points, cond=',1pe12.2)

  write ( *, * ) '#    i       x               numer          anal', &
    '          err          ux (numer and anal)          d/dxerr'

  npmx = 201
  xfun = 0.0
  dxfun = 0.0
  xxf = 0.0
  xdx = 0.0
  do k = 1, npmx
    z = x1 + ( x2 - x1 ) * real ( k - 1 ) / real ( npmx - 1 )
    sol = 0.0
    dsdx = 0.0
    do j = 1, n
      sol = sol + y(j) * g(j,z)
      dsdx = dsdx + y(j) * gx(j,z)
    end do

    anal = fun(z)
    err = anal - sol

    adfdx = dfdx(z)
    errdf = adfdx - dsdx

    if ( abs ( err ) > abs ( xfun ) ) then
      xfun = err
      xxf = z
    end if

    if ( abs ( errdf ) > abs ( dxfun ) ) then
      dxfun = errdf
      xdx = z
    end if
!
!  write (*,3100)k,z,sol,anal,err,dsdx,adfdx,errdf
  write ( io, 3100 ) k,z,sol,anal,err,dsdx,adfdx,errdf
3100  format(i5,1p7e15.5)
  end do

  rewind io

  write (*,3150)xfun,xxf
3150  format(' max. func. error ',1pe15.5,' at ',1pe15.5)
  write ( *,3151)dxfun,xdx
3151  format (' max. df/dx error ',1pe15.5,' at ',1pe15.5)
!
!  write solution at the nodes
!
  open ( io2, file = 'mqle2.dat', status = 'unknown' )
  rewind io2

  write (io2,3060)n,cond
3060  format('# solution ',i5,' points, cond=',1pe12.2)

  write ( *, * ) ' max f,fx error=',xfun, xxf
  write (io2  ,3061)
3061  format('#',4x,'i',7x,'x',15x,'numer',10x,'anal',10x,'err',10x, &
     'ux (numer and anal)',10x,'d/dxerr' &
     ,10x,'c[j]',10x,'curv',10x,'rad',12x,'uxx')

  xfun = 0.0
  dxfun = 0.0
  xxf = 0.0
  xdx = 0.0
  do k = 1, n
    z = xk(k)
    sol = 0.0
    dsdx = 0.0
    do j = 1, n
      sol = sol + y(j) * g(j,xk(k))
      dsdx = dsdx + y(j) * gx(j,xk(k))
    end do

    anal = fun(xk(k))
    err = anal - sol

    adfdx = dfdx(xk(k))
    errdf = adfdx - dsdx

    if ( abs ( err ) > abs ( xfun ) ) then
      xfun = err
      xxf = xk(k)
    end if

    if ( abs ( errdf ) > abs ( dxfun ) ) then
      dxfun = errdf
      xdx = xk(k)
    end if

    ck = abs ( c(k) / y(k) )
    uxx = d2fdx2 ( xk(k) )
    rcur = uxx / ( sqrt ( 1.0 + adfdx**2 ) )**3
    radc = 300.0

    if ( abs ( rcur ) > 1.0 / radc ) then
      radc = 1.0 / abs ( rcur )
    end if

    write (io2 ,3250)k,xk(k),sol,anal,err,dsdx,adfdx,errdf,ck,rcur,radc,uxx
3250  format(i5,1p11e15.5)

  end do

  rewind io2
  close(io2)
!
!   test points
!
  20  continue

  write ( *, * ) ' '
  write ( *, * ) 'Next command?'
  write ( *, * ) ' '
  write ( *, * ) '  Q = quit'
  write ( *, * ) '  G = gnuplot'
  write ( *, * ) '  S = save to a file'
  write ( *, * ) '  X = enter a point X.'

  read ( *, '(a)' ) par

  if ( par(1:1) == 'Q' ) then
    stop
  else if ( par(1:1) == 'S' )then
    write ( *, * ) ' '
    write ( *, * ) 'Enter file name:'
    read ( *, '(a)' ) par
    comnd='cp -p mqle2.dat '// par
    write(*,'(a,a)')' comnd = ', comnd
    call system ( comnd )
    go to 20
  else if ( par(1:1) == 'G' )then
    call gnuplot
    go to 20
  else if ( par(1:1) == 'X' ) then
    read ( par, * ) z
  else
    stop
  end if
!
!  Find the numerical solution at the point Z.
!
   sol = 0.0
   dsdx = 0.0
   do j = 1, n
     sol = sol + y(j) * g(j,z)
     dsdx = dsdx + y(j) * gx(j,z)
   end do

  anal = fun(z)
  err = anal - sol
  adfdx = dfdx(z)
  errdf = adfdx - dsdx

  write ( *,3800)