Discussion :

[Bob_Tom]

bonjour,

j'ai une problème lors de l'écriture d'un programme, je souhaite réalisé x^a, cela marche pour x positif mais ne marche pas pour x négatif sortie 'NaN'?

Ligne 30 du programme je sort x et x^expo et pour les valeur négative de x, x^expo est 'NaN'

Toute mes real sont déclaré en double précision, et j utilise le compilateur gfortran.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
 
Program grounding_line_dinamic
c resouldre le model B de la publication de Scholz 2007
	implicit none
 
	include 'parameter.data'
	include 'physic_parameter.data'
 
	double precision hice(N1),bed_shape(N1),x_gl
	double precision fbed,fbedl
	double precision cst_16,solvxg(N1),ii,solxg(4),tmp,expo,tmp2
	integer i
 
c#############trouver la position de la grounding line#############
c	résolution de l'équation 16
	cst_16=(((A*(rhoi*g)**(n+1.)*(1.-rhoi/rhow)**n)/(C*4.**n))
     #	**(1./(m+1.)))
 
 
	print*,'constante : ',cst_16
	expo=(n+m+3.)/(m+1.)
	print*,expo
	open(1,file='fonction16.data')
	do i=1,N1
	ii=float(i-1)
c	print*,i,ii*dsig*xmax,fbed(ii*dsig*xmax)
 
	tmp=fbedl(ii*dsig*xmax)
 
	tmp2=tmp**expo
	print*,tmp,tmp2
	solvxg(i)=cst_16*(rhow/rhoi)**expo
     #	*tmp2-acu*ii*dsig*xmax
 
	write(1,*)ii*dsig*xmax,solvxg(i)
	enddo
	close(1)
 
c	call rootvect(solvxg,solxg)
 
	do i=1,4
	print*,solxg(i)
	enddo
c	fin du programme
	end
c--------------------------------------------------------------------
	function fbedl(x)
	implicit none
 
	double precision fbedl,x
 
	fbedl=-729.d0+1.11d0*x
 
 
	return
	end function

pour les includes

parameter.data

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
 
c	N1 step decompositon first part
 
	integer N1
	parameter(N1=100)
	double precision dsig,xmax
 
	parameter(xmax=1450.,dsig=1./float(N1-1))

physic_parameter.data

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
	double precision m,n,rhoi,rhow,g,C,A,acu
 
	parameter(g=9.8,C=7.624*10**6,rhoi=900.,rhow=1000.)
	parameter(acu=0.3,A=4.23*10.**(-25.))
	parameter(n=3.,m=1./n)

quand je fait un petit progamme comme celui-ci pas de problème

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
	program test
	implicit none
	double precision B
 
	print*,-729.123**(4.75)
	end

Merci de votre aide

[Sylvain Bergeron]

Ça ne marche pas car c'est ce que dit la norme. Un réel négatif à un exposant réel est illégal. Certains compilateurs acceptent tout de même d'effectuer le calcul si l'exposant réel correspond à un entier (ex: 2.0).

Pourquoi ? Tu essayeras de tracer la courbe -2^x pour x < 0...

[Bob_Tom]

je ne suis pas sur que cela vienne de ca car dans le petit programme cela ne pose pas de soucie
de plus ici x^a ou 'a' est real positif

[Marlan]

Salut !

Pourquoi ne pas rajouter une condition logique du type :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
 
if(a<0d0) then
CALC=1d0/(x**(-a))
else
CALC=x**a
end if

?

[Bob_Tom]

Le problème n'est pas a négatif mais pour x négatif ! x^a dans mon problème 'a' a une valeur fixe 29/4 mais que je définie comme un réel et pas un rationelle.

j'ai plus ou moins contournez le problème avec

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
	tmp=fbed(ii*dsig*xmax)
	if (tmp.lt.(0.)) then
	tmp2=-1.*(-1.*tmp)**expo	
	else
	tmp2=tmp**expo
	endif

mais ce que je ne comprend pas c'est pourquoi cela pose un problème ici alors que dans le mini prog que j'ai fait il n'y en a pas (cf premier post)

[Marlan]

Re-salut,

Ok, j'avais mal compris la question.

Le problème ici est que tu cherches à calculer avec un exposant non entier --> une racine (tordue certes, mais une racine quand même). Or, une racine ne peut jamais être négative dans l'espace réel.

En somme, si tu écris (-X)**(29/4), ou bien X**(29/4) avec X<0, tu calcules une racine de nature complexe et non pas réelle.

De fait, dans ton code, il va être nécessaire d'isoler le signe de X via quelque chose comme :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
 
SGN=X/abs(X)
CALC=SGN*abs(X)**a

Et là ça devrait marcher.

En espérant avoir répondu à ta question cette fois ci,

Marlan