I am trying to illustrate how to pass a function to the Newton Raphson procedure. I manage to execute a very simple function (called unefonctionsee below), but it does not work with a function that has parameters. This second highlight is called gaussienne, and it takes one argument, x and two optional arguments, muand sig. I called her so in my procedure Raphson newton: f(x). What is strange to me is that at run time the program does as if the optional parameters were present sigand mu, but they do not ... So, I donโt understand ...
Here is a module that contains functions
module fonction
implicit none
double precision :: f_sigma = 1.d0, f_mu = 0.d0
double precision, parameter :: pi = 3.14159265359d0
contains
double precision function unefonction(x)
implicit none
double precision, intent(in) :: x
unefonction = (exp(x) - 10.) / (x + 2.)
end function unefonction
double precision function gaussienne(x, mu, sig)
implicit none
double precision, intent(in) :: x
double precision, intent(in), optional :: mu, sig
double precision :: norme, moy, sigma
if (present(sig)) then
write(*,*)"sig present"
sigma = sig
else
sigma = f_sigma
end if
if (present(mu)) then
write(*,*)"mu present"
moy = mu
else
moy = f_mu
end if
norme = 1.d0 / (sigma * sqrt(2.d0 * pi))
gaussienne = norme * exp(-(x - moy)**2 / (2.d0 * sigma**2))
end function gaussienne
end module fonction
Here is the module that contains the newton raphson procedure
module rechercheRacine
implicit none
contains
subroutine newtonRaphson(racine, f, eps, cible)
implicit none
double precision, intent(inout) :: racine
double precision, intent(in), optional :: cible, eps
double precision, external :: f
integer :: compteur
double precision :: xold, xnew, delta, valcible
double precision :: threshold, fprim, fdex
if (present(eps)) then
threshold = eps
else
threshold = 1.d-10
end if
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(*,*) "--------------------------------------------------------"
compteur = 0
delta = 1.d0
xold = racine
write(*, '(i4,4e16.6)') compteur, f(xold), xold, 0., threshold
do while (delta > threshold .and. compteur <= 100)
fdex = f(xold) - valcible
fprim = (f(xold + threshold) - f(xold - threshold)) / (2.d0 * threshold)
xnew = xold - fdex / fprim
delta = abs(xnew - xold)
compteur = compteur + 1
write(*, '(i4,4e16.6)') compteur, fdex, xnew, delta, threshold
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
use rechercheRacine
use fonction
implicit none
double precision :: racine, eps, cible
racine = 1.d0
call newtonRaphson(racine, unefonction)
racine = 1.d0
eps = 1.d-14
cible = 10.d0
call newtonRaphson(racine, unefonction, eps, cible)
f_sigma = 2.d0
f_mu = 5.d0
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 , , sig mu, , .