Discussion :

[herkitz]

Bonjour,

Je souhaite calculer l'intégrale (sur un triangle) suivante :

où la fonction phi est définie par :

J'utilise pour ce faire le code suivant :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
 
Phi=function(s,d){
	Phi= (exp(d)-exp(s))/(exp(d)-1)
	return(Phi)
}
Integrande2=function(s,d){
	 Res=max(0,(Phi(s,d)-1/s^5)^2)
	 return(Res)
}
INtegrale=function(){ integrate(function(y) { sapply(y, function(y) { integrate(function(x) Integrande2(x,y), 0, y)$value   }) }, 0, 200)}
INtegrale()

J'ai cependant cette erreur:

Code :

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

1
2
3
 
Erreur dans integrate(function(x) Integrande2(x, y), 0, y) : 
  evaluation of function gave a result of wrong length

Auriez vous une idée?

Merci d'avance pour votre aide.

[dev_ggy]

Bonjour,

Peux-tu indiquer le package que tu utilises ?

Cordialement.

[juliatheric]

Bonjour,
Dans INtegrale, à l'intérieur de sapply, dans le premier argument de integrate, la longueur du vecteur en sortie (le résultat de l'application de Integrande2 sur x et y) doit être égale à la longueur du vecteur en entrée (x). C'est ce qui est spécifié dans la doc de integrate :

f : an R function taking a numeric first argument and returning a numeric vector of the same length. Returning a non-finite element will generate an error.

Une fois ceci corrigé, il y aura un autre problème : l'intégrale serait infinie. Ce pourrait être dû à la manière dont les bornes sont définies dans Integrande2. Alors le code incluant toutes les corrections :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
> Phi <- function(s,d) (exp(d)-exp(s)) / (exp(d)-1)
> h <- 0
> Integrande2 <- function(h, s) max(h, (Phi(s,d) - (1/(s^5)))^2)
> INtegrale <- function(){
	integrate(function(y) {
		sapply(y, function(y)
			integrate(function(x) {
				w <- Integrande2(x,y)
				rep(w, length(x))
				},0, y)$value)
		}, h, 200)
}
 
> INtegrale()

2.396808e+171 with absolute error < 2.4e+166

À voir si cela correspond au résultat attendu.