Discussion :

[Hallu33]

Bonjour, j'ai un petit souci. Je sais intégrer une fonction f(r) avec quad(f(r), a, b) a et b étant les bornes. Mon souci est que j'ai une fonction f(r) qui est en fait de la forme :

f(r)=g(r)*int(h(r,s), s = 0..1)

Il faudrait donc que je calcule une intégrale par rapport à "s" et dont le résultat soit fonction de "r", puis ensuite que j'intègre le tout par rapport à "r" (impossible de trouver une expression exacte, c'est le produit d'une fonction de Bessel et d'un cos). Mais quad ne peut me retourner autre chose qu'un nombre, et ne peut me retourner une fonction dépendant de "r".

Le seul moyen que je vois c'est de créer un vecteur r=[0:1e-3:1] par exemple, et calculer int(h(r,s), s= 0..1) pour chaque valeur de r. Mais dans ce cas comment faire à nouveau quad de quelque chose pour calculer l'intégrale de f(r) qui ne sera alors pas une fonction mais un vecteur de longueur 1000 dans mon cas, sachant que je n'ai évidemment pas envie d'interpoler ces 1000 points pour créer une fonction continue intégrable ? Ce que je recherche c'est une façon automatique de relier ces 1000 points entre eux (faire 999 segments continus mis bout à bout) et calculer l'aire dessous cette nouvelle courbe continue. Y a-t-il une option de quad ou de tout autre système de calcul d'intégrale matlab qui fait ça automatiquement (on lui rentre un vecteur, il interpole tout seul et fait le calcul de l'intégrale) ?

Merci d'avance de votre aide.

[grem1]

bonjour,

dblquad peut être

[FR119492]

Salut!

Ton problème semble provenir du fait que le mot intégrale a deux significations distinctes:

• L'intégrale indéfinie, appelée aussi primitive d'une fonction f(x) est une fonction g(x) dont f(x) est la dérivée; si tu connait l'expression de f(x), la recherche de l'expression analytique de g(x) est du calcul formel, que tu peux faire par exemple avec Maple; tu peux ensuite la tabuler si tu veux la représenter graphiquement.
• L'intégrale définie d'une fonction est une valeur unique qui correspond à la surface au-dessous du graphe de f(x); sa détermination est du calcul numérique.

Ces deux notions sont évidemment liées, mais ce n'est pas une raison pour les confondre.

A part ça, ce n'est pas une question de Matlab mais de mathématiques.

Jean-Marc Blanc

[Hallu33]

Justement calculer la primitive de h(r,s) ou de f(r) est impossible analytiquement, c'est ce que j'ai précisé clairement dans mon sujet. Voilà pourquoi je dois utiliser quad pour les deux.

[FR119492]

Salut!
Est-ce que ton problème n'est pas simplement de calculer une intégrale double?
Jean-Marc Blanc

[Hallu33]

ben je dois en fait calculer un truc du type :

int(g(r)*int(h(r,s), s=0..a), r=0..inf)

Si h était une fonction dont on pouvait calculer la primitive de manière exacte, aucun problème mais c'est impossible. Par conséquent, on ne peut utiliser un calcul numérique de int(h(r,s), s=0..a) autrement qu'à r fixé pour une valeur donnée. L'idée est donc de calculer g(r)*int(h(r,s), s=0..a) pour r allant de 0 à une grande valeur avec un pas suffisamment fin pour pouvoir ensuite créer vraiment la fonction continue f(r)=g(r)*int(h(r,s), s=0..a) par une interpolation des valeurs pour chaque r calculées précédemment, puis d'intégrer f(r) de 0 à l'infini. Ma question je le répète est donc y a-t-il moyen d'automatiser le calcul de l'intégrale de f(r) sans devoir faire une interpolation à la main. Du genre on rentre le vecteur r, la fonction f, et ça interpole tout seul et calcule l'intégrale numérique de la fonction interpolée tout seul.

[FR119492]

Salut!

int(g(r)*int(h(r,s), s=0..a), r=0..inf)

Si tu nous disais d'où vient cette formule, on pourrait peut-être mieux t'aider: il vaut souvent mieux commencer par le commencement.
Jean-Marc Blanc

[Hallu33]

Il s'agit du champ acoustique rayonné par un transducteur focalisé d'ouverture numérique 1.
Le champ est calculé dans l'approximation de Huygens-Fresnel, sous la forme d'un potentiel vitesse qui dépend de r et de z. Je dois dériver ce potentiel par rapport à z, et intégrer le produit de ce potentiel et de sa dérivée par rapport à r.

Voici le champ Psi(r,z) en pièce jointe que je dois ensuite intégrer par rapport à r.

[salseropom]

ben je dois en fait calculer un truc du type :

int(g(r)*int(h(r,s), s=0..a), r=0..inf)

Si h était une fonction dont on pouvait calculer la primitive de manière exacte, aucun problème mais c'est impossible. Par conséquent, on ne peut utiliser un calcul numérique de int(h(r,s), s=0..a) autrement qu'à r fixé pour une valeur donnée. L'idée est donc de calculer g(r)*int(h(r,s), s=0..a) pour r allant de 0 à une grande valeur avec un pas suffisamment fin pour pouvoir ensuite créer vraiment la fonction continue f(r)=g(r)*int(h(r,s), s=0..a) par une interpolation des valeurs pour chaque r calculées précédemment, puis d'intégrer f(r) de 0 à l'infini. Ma question je le répète est donc y a-t-il moyen d'automatiser le calcul de l'intégrale de f(r) sans devoir faire une interpolation à la main. Du genre on rentre le vecteur r, la fonction f, et ça interpole tout seul et calcule l'intégrale numérique de la fonction interpolée tout seul.

Salut, tu viens de dire l'algo à utiliser :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
 
 
step = 1e-3;
rmax = 1000;
iter = 0;
vectr = 0 : step : rmax;
res = zeros(size(vectr));
 
for r = vectr
  1) calculer aux = int(h(r,s), s=0..a) --> voir la méthode de Simpson
  2) iter = iter + 1;
  3) res(iter) = g(r) * aux;
end % for r = vectr

Pour calculer l'intégrale de f(r), comme tu ne connais ta fonction f que sur qq points (fonction discrétiser), tu peux utiliser la méthode des trapèzes qui est exacte pour les polynômes de degré 1 (interpolation linéaire en somme). En matlab tu utilises la fonction trapz. Ca donne :

Code :

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

1
2
 
resultat = trapz(vectr,res);

Je n'ai pas testé ce code, mais ça devrait marcher.

[FR119492]

Salut!

sous la forme d'un potentiel vitesse

Petite incohérence: un potentiel est par nature un champ scalaire et une vitesse est par nature un champ vectoriel.
Jean-Marc Blanc

[Hallu33]

En physique on parle tous les jours de potentiel vitesse (velocity potential en anglais). C'est juste le potentiel dont découle le vecteur vitesse.

Sinon merci pour le code j'ignorais l'existence de cette fonction trapèze. Enfin finalement je suis tout de même parti sur une approximation gaussienne du champ bessélien, le calcul de la double intégrale est trop long en terme de temps.

[FR119492]

potentiel vitesse

En français, le potentiel du champ de vitesse, car une vitesse s'exprime en mètres par seconde.
Ce potentiel obéit à une équation de d'Alembert avec des conditions aux limites données. Il suffit d'intégrer cette équation après avoir choisi le système de coordonnées le mieux adapté.
Jean-Marc Blanc

[Loïc B.]

Je repars du premier post, parce que dans ce cas, ton énoncé était clair, nul besoin de savoir si t'intègre des potentiels, des intégrales définies, indéfinies, ou des patates !

Tu veux intégrer f(r) :

int(f(r),r) = int(g(r)*int(h(r,s),s) ,r)

1 - Mathématique :
ta fonction g(r) est une constante par rapport à la variable s, puisqu'indépendante de celle-ci. Tu peux donc la rentrer dans la seconde intégrale :

int(f(r),r) = int(int(g(r)dr* h(r,s)ds,s),r)

Et là magie !! Une intégrale double.

2 - Matlab
tu as donc une solution toute trouvée (et soufflée par grem1)

dblquad(@(r,s)g(r).*h(r,s),rmin,rmax,smin,smax)

Attention, on n'oublie pas le point dans la multiplication.

Et voilà.

PS : je connaissais pas dblquad.
PPS : je me réserve le droit de me tromper.

[Hallu33]

Merci ça m'aide pour résoudre une partie du problème. La suite est plus compliquée cependant : en fait f dépend aussi de z, soit f(r,z). Je dérive f par rapport à z. Ca me donne en gros du g(r)*[int(h(r,s),s)^2+int(h(r,s),s)*int(i(r,s),s)]

Comment calculer du coup l'intégrale sur r de g(r)*int(h(r,s),s)*int(i(r,s),s) ?

[edit] autant pour moi je crois avoir trouvé, il suffit de faire une "fausse" intégrale triple, en lui disant que i(r,s) dépend en fait de i(r,s') mais avec les mêmes bornes.

[Hallu33]

Bon ça marche, problème : le calcul prend une minute environ, et je dois faire le calcul pour 9 triples intégrales, le tout sur 5000 cycles, soit 45 000 calculs =(

[nahouto]

Bon ça marche, problème : le calcul prend une minute environ, et je dois faire le calcul pour 9 triples intégrales, le tout sur 5000 cycles, soit 45 000 calculs =(

tu n'as qu'à utiliser une grille de calcul (je crois que le CEA de Saclay en prête/loue)

[Hallu33]

lol merci mais je doute que le CEA prête ses grilles de calcul sans faire sa petite enquête sur le calcul en question, j'ai un pote qui y a fait un post-doc c'était l'horreur pour y entrer, avec des questions aussi débiles que "avez-vous déjà été activiste syndicaliste ?".

[nahouto]

lol merci mais je doute que le CEA prête ses grilles de calcul sans faire sa petite enquête sur le calcul en question, j'ai un pote qui y a fait un post-doc c'était l'horreur pour y entrer, avec des questions aussi débiles que "avez-vous déjà été activiste syndicaliste ?".

Si tu as une connaissance qui y bosse, tu peux tenter.
Parce que sinon, tu devras attendre un bon bout de temps avant de voir tes résultats