2013-02-25 103 views
2

我嘗試說明如何將函數傳遞給Newton Raphson過程。我用一個非常簡單的函數(稱爲unefonction見下文)成功,但它不適用於具有參數的函數。這第二個功能稱爲gaussienne,它有一個參數x和兩個可選參數musig。在我的牛頓raphson程序中,我以這種方式稱作功能:f(x)。對我來說奇怪的是,在執行過程中,程序的作用好像有可選參數sigmu,但它們不會...因此我不明白...將函數作爲參數傳遞給子例程時出現分段錯誤

這裏是包含模塊功能

module fonction 

    implicit none 

    ! parametre pour la gaussienne 
    double precision :: f_sigma = 1.d0, f_mu = 0.d0 

    ! pi accessible uniquement en interne 
    double precision, parameter :: pi = 3.14159265359d0 

    contains 

    double precision function unefonction(x) 
     ! fonction : unefonction 
     ! renvoie 
     ! $\frac{e^x - 10}{x + 2}$ 

     implicit none 

     ! arguments 
     double precision, intent(in) :: x 

     unefonction = (exp(x) - 10.)/(x + 2.) 

    end function unefonction 

! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 

    double precision function gaussienne(x, mu, sig) 
     ! fonction gaussienne 
     ! utilise les parametres definis dans le module si 
     ! mu et sig ne sont pas passes en argument 

     implicit none 

     ! arguments 
     double precision, intent(in)   :: x 
     double precision, intent(in), optional :: mu, sig 

     ! variables locales 
     double precision :: norme, moy, sigma 

     ! sigma 
     if (present(sig)) then 
      write(*,*)"sig present" 
      sigma = sig 
     else 
      sigma = f_sigma 
     end if 

     ! mu 
     if (present(mu)) then 
      write(*,*)"mu present" 
      moy = mu 
     else 
      moy = f_mu 
     end if 

     ! calcul de la gaussienne 
     norme = 1.d0/(sigma * sqrt(2.d0 * pi)) 
     gaussienne = norme * exp(-(x - moy)**2/(2.d0 * sigma**2)) 

    end function gaussienne 

end module fonction 

這裏是包含牛頓拉夫遜過程

module rechercheRacine 

    implicit none 

    contains 

     subroutine newtonRaphson(racine, f, eps, cible) 

      ! recherche l'antecedant de cible 

      implicit none 

      ! arguments 
      double precision, intent(inout)  :: racine 
      double precision, intent(in), optional :: cible, eps 

      ! fonction dont on cherche la racine 
      double precision, external :: f 

      ! variables locales 
      integer   :: compteur 
      double precision :: xold, xnew, delta, valcible 
      double precision :: threshold, fprim, fdex 

      ! precision 
      if (present(eps)) then 
       threshold = eps 
      else 
       threshold = 1.d-10 
      end if 

      ! valeur cible 
      if (present(cible)) then 
       valcible = cible 
      else 
       valcible = 0.d0 
      end if 

      write(*,*) "--------------------------------------------------------" 
      write(*,*) "      NEWTON RAPHSON" 
      write(*,*) "--------------------------------------------------------" 
      write(*,"('x0 = ',e16.6)") racine 
      write(*,"('seuil = ',e16.6)") threshold 
      write(*,"('cible = ',e16.6)") valcible 
      write(*,*) "--------------------------------------------------------" 
      write(*,*) "      ITERATIONS" 
      write(*,*) "--------------------------------------------------------" 

      ! initialisation 
      compteur = 0 
      delta = 1.d0 
      xold  = racine 
      write(*, '(i4,4e16.6)') compteur, f(xold), xold, 0., threshold 

      ! iterations 
      do while (delta > threshold .and. compteur <= 100) 

       ! calcul de la fonction en xold 
       fdex = f(xold) - valcible 

       ! calcul de la derivee numerique 
       fprim = (f(xold + threshold) - f(xold - threshold))/(2.d0 * threshold) 

       ! application de l'iteration de Newton Raphson 
       xnew = xold - fdex/fprim 
       delta = abs(xnew - xold) 
       compteur = compteur + 1 

       ! affichage de la convergence 
       write(*, '(i4,4e16.6)') compteur, fdex, xnew, delta, threshold 

       ! mise a jour de xstart 
       xold = xnew 
      end do 

      if (delta < threshold) then 
       racine = xnew 
       write(*, *) '--------------------------------------------------------' 
       write(*, *) '      CONVERGE' 
       write(*, *) '--------------------------------------------------------' 
       write(*, *) 'A la convergence demandee, une solution est:' 
       write(*, "('x = ',e20.10,' f(x) = ', e20.10)") racine, f(racine) 
       write(*, *) 
      else 
       write(*, *) '--------------------------------------------------------' 
       write(*, *) '      NON CONVERGE' 
       write(*, *) '--------------------------------------------------------' 
      end if 

     end subroutine newtonRaphson 

end module rechercheRacine 

這裏是主程序中的模塊:

program main 

    ! contient la subroutine newtonRaphson 
    use rechercheRacine 

    ! contient la fonction 
    use fonction 

    implicit none 

    double precision :: racine, eps, cible 

    ! appel de la subroutine newtonRaphson 
    ! sans la valeur cible : cible (defaut = 0) 
    ! sans la precision : eps (defaut 1d-10) 
    racine = 1.d0 
    call newtonRaphson(racine, unefonction) 

    ! -------------------------------------------------------- 

    ! appel de la subroutine newtonRaphson 
    ! avec pour cible 10 
    racine = 1.d0 
    eps = 1.d-14 
    cible = 10.d0 
    call newtonRaphson(racine, unefonction, eps, cible) 

    ! -------------------------------------------------------- 

    ! parametre de la gaussienne 
    f_sigma = 2.d0 
    f_mu = 5.d0 
    ! appel de la subroutine newtonRaphson 
    ! passage des arguments sous la forme clef = valeur 
    cible = 0.1d0 
    racine = 2.d0 
    call newtonRaphson(cible = cible, f = gaussienne, racine = racine) 

end program main 

主程序適用於稱爲unefonction的功能,但它不適用於gaussienne功能。

以下是錯誤消息:

Program received signal SIGSEGV: Segmentation fault - invalid memory reference. 

Backtrace for this error: 
#0 0x7F1B6F5890F7 
#1 0x7F1B6F5896D4 
#2 0x7F1B6EEEB49F 
#3 0x4009D2 in __fonction_MOD_gaussienne at mod_fonction.f90:54 
#4 0x40104D in __rechercheracine_MOD_newtonraphson at mod_racine.f90:59 
#5 0x4016BA in MAIN__ at main.f90:40 
Erreur de segmentation (core dumped) 

我認爲invalid memory reference是由於這樣的事實,該方案確實好像可選參數sigmu存在,從而尋找他們,而他們不是。

回答

6

是的,問題是你傳遞的函數確實需要三個參數而不是一個。如果您在子程序newtonRaphson

double precision, external :: f 

改變external聲明f以顯式接口(描述,如何真正使用它):

interface 
    double precision function f(x) 
    double precision, intent(in) :: x 
    end function f 
end interface 

你的代碼甚至不會編譯由於參數數量不匹配。

他們不同的方式來傳遞「參數」的功能f這是從常規newtonRaphson叫:

  • 你可以期望f有兩個intent(in)參數,而不是一個:附加到值x ,它也可能需要一個真實的數組,它可以是任意大小的,並且可以包含必要的參數。這需要以下接口:

    interface 
        double precision function f(x, params) 
        double precision, intent(in) :: x 
        double precision, intent(in) :: params(:) 
        end function f 
    end interface 
    

    這些功能,它不需要參數(如unefonction)將只是沒有使用第二個參數的內容,而其他(如gaussienne)將持有它們的參數從它。

  • 您可以使newtonRaphson預期給定的可擴展類型(class)和一個類型綁定過程返回給定x值的值。然後,您可以創建這種類型的abritrary擴展,它可以基於存儲爲擴展類型中的字段的一些參數來計算給定x值的值。然後該程序可能看起來像下面(我剝幾個部分),但它需要的Fortran編譯器2003:

    module rechercheRacine 
        implicit none 
    
        type, abstract :: calculator 
        contains 
        procedure(getvalue_iface), deferred :: getvalue 
        end type calculator 
    
        interface 
        double precision function getvalue_iface(self, x) 
         import :: calculator 
         class(calculator), intent(in) :: self 
         double precision, intent(in) :: x 
        end function getvalue_iface 
        end interface 
    
    contains 
    
        subroutine newtonRaphson(racine, f, eps, cible) 
        double precision, intent(inout)  :: racine 
        class(calculator), intent(in)   :: f 
        double precision, intent(in), optional :: cible, eps 
    
        do while (delta > threshold .and. compteur <= 100) 
         fdex = f%getvalue(xold) - valcible 
         : 
        end do 
    
        end subroutine newtonRaphson 
    
    end module rechercheRacine 
    
    
    module fonction 
        use rechercheRacine 
        implicit none 
    
        type, extends(calculator) :: unefonction 
        contains 
        procedure :: getvalue => unefonction_getvalue 
        end type unefonction 
    
        type, extends(calculator) :: gaussienne 
        double precision :: mu = 0.0d0, sigma = 1.0d0 
        contains 
        procedure :: getvalue => gaussienne_getvalue 
        end type gaussienne 
    
    contains 
    
        double precision function unefonction_getvalue(self, x) 
        class(unefonction), intent(in) :: self 
        double precision, intent(in) :: x 
        unefonction_getvalue = (exp(x) - 10.)/(x + 2.) 
        end function unefonction_getvalue 
    
        double precision function gaussienne_getvalue(self, x) 
        class(gaussienne), intent(in) :: self 
        double precision, intent(in) :: x 
    
        : 
        gaussienne_getvalue = norme * exp(-(x - moy)**2/(2.d0 * self%sigma**2)) 
    
        end function gaussienne_getvalue 
    
    end module fonction 
    
    
    program main 
        use rechercheRacine 
        use fonction 
        implicit none 
    
        type(unefonction) :: fone 
        type(gaussienne) :: fgauss 
    
        : 
        call newtonRaphson(racine, fone) 
        call newtonRaphson(cible = cible, f = fgauss, racine = racine) 
    
    end program main 
    
+0

尼斯的答案。實際上,外部聲明只是我原諒刪除的一個測試。我認爲如果一個函數/子程序在一個模塊中,則會自動創建一個顯式接口。我錯了嗎 ?但我明白,外部聲明將刪除接口。 我認爲參數sig和mu的可選聲明是我沒有在編譯時出錯的原因。 – Ger 2013-02-25 22:11:51

相關問題