Discussion :

[Beltharion]

Boujour à toutes et à tous,

Le programme que je suis en train d'écrire à pour but de calculer un champ de vitesse dans un système eu peu compliquer. J'ai résolu ce calcul de façon analytique en terme de transformée de Hankel.
Il me faut donc inverser le calcul analytique en utilisant une procédure numérique. Et c'est là le problème !!!

Avant de parler du programme à proprement parler, il me faut d'abord parler de la résolution analytique.

Comme je l'ai dit précédemment, le calcul est résolu dans l'espace de Hankel. Le-dit calcul consiste à résoudre un système différentiel avec certaines conditions aux limites. Il y a donc deux étapes importantes dans le calcul :
1) le calcul des paramètres de résolutions du système différentiel. Il s'agit dans mon cas d'un système linéaire de 6 équations à 6 inconnues. Normalement, cela ne devrait pas poser de problèmes.
2) Le calcul du champ de vitesse dans l'espace de Hankel avec les paramètres issus de 1)

C'est là qu'intervient la résolution numérique ...
Donc le programme pourrait se structurer en 3 parties :
1) Résolution des conditions limites
2) Injection du résultat de 1) pour calculer le champ de vitesse dans l'espace Hankel
3) Inversion de la transformée de Hankel pour obtenir le champ de vitesse dans l'espace des réels

Cependant en faisant une transformée une transformée de Hankel on associe à la variable réelle (ici r) une variable de dimension inverse (ici k=nombre d'ondes) par un calcul intégrale évalué entre 0 et +inf. Pour repasser dans l'espace réel, on inverse la transformation toujours entre 0 et +inf. C'est là que réside mon problème.
En effet la matrice qui sort de la résolution des conditions limites automatiquement dépendante de k. Donc pour résumer j'ai :
M(k)*X=N(k)

Comme vous pourrez le constater, j'ai tenté de fonctionnaliser les matrices M, N et X, vu qu'on ne peut pas calculer une boucle infinie. De cette façon j'espérais pouvoir passer de l'étape 1) à 3) en utilisant la fonction quad(fct,0,Inf) pour calculer brutalement le transformée inverse avec le noyau de Kernel associé (partie III du programme). Cependant, je comprend bien que la notation X(i,j) ne fonctionne pas étant donné que X est fonction de k.

Je n'ai absolument aucune idée de comment résoudre ce problème à cause du caractère impropre de la transformée.

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
 
%% Parameters
 
T0=20;
delta_T0=10;
Tmax=T0+delta_T0;
sigma_vs_T=1E-4;
wth=5E-5; 
nb_pts=1000;
if strcmp(mode_exe,'yes')
    r=linspace(0,1E-3,nb_pts);
    z=linspace(0,7E-4,nb_pts);
else
    r=linspace(0,1E-3,nb_pts);
    z=linspace(0,7E-4,nb_pts);
end
 
%% Calculation
 
Tk=@(k) Tmax*((wth^2)/2).*exp(-((wth.*k).^2)/4);
for i=1:length(r)
    integrand=@(k)(-k*Tmax*(wth^2)/2)*exp(-((wth*k)^2)/4)*k*besselj(1,k*r(i));
    gradTr(i)=quad(integrand,0,inf);
end
 
M=@(k)[-2 -1/(k*H1) -1/(k*H1) 2 1/(k*H2) 1/(k*H2) \
(exp(-k*H1)-exp(k*H1)) -exp(-k*H1) -exp(k*H1) 0 0 0 \
0 0 0 (exp(k*H2)-exp(-k*H2)) exp(k*H2) exp(-k*H2)\
(-exp(-k*H1)-exp(k*H1)) -((1-k*H1)/(k*H1))*exp(-k*H1) -((1+k*H1)/(k*H1))*exp(k*H1) 0 0 0 \
0 0 0 (-exp(k*H2)-exp(-k*H2)) -((1+k*H2)/(k*H2))*exp(k*H2) -((1-k*H2)/(k*H2))*exp(-k*H2) \
0 -2*(eta1/H1) 2*(eta1/H1) 0 2*(eta2/H2) -2*(eta2/H2)];
%%% Matrix N
N=@(k) [0;0;0;0;0;-sigma_vs_T*k*Tmax*((wth^2)/2)*exp(-((wth*k)^2)/4)];
%%% Calculation X=inv(M).N
X=@(k) inv(M)*N;
 
%% %%%% Partie III : Velocity field calculation : inverse hankel transform
%% Calculation of urr
for i=1:length(r)
  for j=1:length(z)
       urr_1(i,j)=quad((-X*(1,1)*(exp(k*z(j))+exp(-k*z(j)))-(X(2,1)/(k*H1))*(1+k*z(j))*exp(k*z(j))+(X(3,1)/(k*H1))*(1-k*z(j))*exp(-k*z(j)))*k*besselj(1,k*r(i)),0,Inf);
 
       urr_2(i,j)=quad((-X(4,1)*(exp(k*z(j))+exp(-k*z(j)))-(X(5,1)/(k*H2))*(1+k*z(j))*exp(k*z(j))+(X(6,1)/(k*H2))*(1-k*z(j))*exp(-k*z(j)))*k*besselj(1,k*r(i)),0,Inf);       
  end
end

Je remercie par avance toute personne qui pourrait m'apporter sons aide. Si d'aventure vous avez des suggestions pour optimiser ce calcul, je suis preneur.


Beltharion

[Invité]

Bonsoir,

Pour rester avec l'utilisation de fonctions anonymes, on pourrait utiliser la fonction de bas niveau subsref comme ceci :

Code :

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

Xkij = @(k, i,j) subsref(X(k),struct('type','()' ,'subs',{{i,j}}));

Puis l'utiliser dans la fonction quad comme ceci :

Code :

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

urr_1(i,j)=quad(@(k)... Xkij(k,2,1) ... ))

Mais cela reviendrait à re-calculer X(k) à chaque apparition de Xkij.

Aussi je te conseillerais plutôt d'utiliser une "vraie" fonction non anonyme qui prendrais la forme suivante, avec M, N, et X qui seraient alors de simples tableaux :

Code :

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

1
2
3
4
5
6
function res = maFonction(k, z, <les constantes>)
 
M = [  ... ];
N = [ ... ];
X = inv(M)*N;
res = (-X*(1,1)*(exp(k*z(j))+exp(-k*z(j)))-(X(2,1)/(k*H1))*(1+k*z(j))*exp(k*z(j))+(X(3,1)/(k*H1))*(1-k*z(j))*exp(-k*z(j)))*k*besselj(1,k*r(i));

Et qui serait utilisée comme suit :

Code :

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

       urr_1(i,j)=quad(@(k) maFonction(k, z, <constantes>)

Avec les <constantes> étant H1, H2, eta1, eta2, sigma_vs_T, Tmax et wth si j'en ai pas oublié.

[Beltharion]

Je te remercie pour ta réponse détaillée,

Vu que le calcul de champ de vitesse est assez long, j'ai tendance à penser que l'utilisation d'une fonction non anonyme me ferai gagner du temps.
Dans ce cas je crée la fonction maFonction(k, z, <les constantes>).

Mais ensuite quand je l'appelle dans la fonction quad de la façon suivante :

Code :

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

urr_1(i,j)=quad(@(k) maFonction(k, z, <constantes>)

Le calcul peut-il être évalué sur l'intervalle de définition de la transformée ?

Code :

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

urr_1(i,j)=quad(@(k) maFonction(k, z, <constantes>,0,Inf)

[Invité]

La fonction quad non, mais la fonction quadgk, qui s'utilise de la même façon, si

[Beltharion]

Je suis en train d'implémenter tout ça avec la fonction quadgk.

Par contre je ne testerai que demain

En tout cas merci pour ta réponse rapide et précise.

[Beltharion]

J'ai implémenté la fonction mais je dois faire face à une autre erreur.

Dans un premier voilà le code de la fonction :

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
57
58
59
 
% Calculation of the parameters from the Boundary conditions for urk_1
 
% This function compute the inverse hankel tranform defined between 0 and +inf
 
% K         variable in Hankel space
% Z         variable for height defined between 0 and H
% R         variable in real space associated to k
% Other parameters are constant defined in parameters.m. These constants are not required to 
% run the function
 
 
function [res]=BC_r1(K,Z,R,H1,H2,eta1,eta2,Tmax,sigma_vs_T,wth)
  if nargin < 1|| isempty(K)
    disp('Not enough arguments')
  end
  if nargin < 2|| isempty(Z)
    disp('Not enough arguments')
  end
  if nargin < 3|| isempty(R)
    disp('Not enough arguments')
  end
  if nargin < 4 || isempty(H1)
    load parameters
    H1=parameters.H1;
    H2=parameters.H2;
    eta1=parameters.eta1;
    eta2=parameters.eta2;
    Tmax=parameters.Tmax;
    sigma_vs_T=parameters.sigma_vs_T;
    wth=parameters.wth;
  end
%  K=1;
%  R=1;
%  Z=1;
%  load parameters
%  H1=parameters.H1;
%  H2=parameters.H2;
%  eta1=parameters.eta1;
%  eta2=parameters.eta2;
%  Tmax=parameters.Tmax;
%  sigma_vs_T=parameters.sigma_vs_T;
%  wth=parameters.wth;
 
  % Matrix definition
  L1=[-2 -1/(K*H1) -1/(K*H1) 2 1/(K*H2) 1/(K*H2)];
  L2=[(exp(-K*H1)-exp(K*H1)) -exp(-K*H1) -exp(K*H1) 0 0 0];
  L3=[0 0 0 (exp(K*H2)-exp(-K*H2)) exp(K*H2) exp(-K*H2)];
  L4=[(-exp(-K*H1)-exp(K*H1)) -((1-K*H1)/(K*H1))*exp(-K*H1) -((1+K*H1)/(K*H1))*exp(K*H1) 0 0 0];
  L5=[0 0 0 (-exp(K*H2)-exp(-K*H2)) -((1+K*H2)/(K*H2))*exp(K*H2) -((1-K*H2)/(K*H2))*exp(-K*H2)];
  L6=[0 -2*(eta1/H1) 2*(eta1/H1) 0 2*(eta2/H2) -2*(eta2/H2)];
  M=[L1;L2;L3;L4;L5;L6];
  N=[0;0;0;0;0;-sigma_vs_T*K*Tmax*((wth^2)/2)*exp(-((wth*K)^2)/4)];
  % Calculation
  [L,U]=lu(M);
  temp=pinv(L)*N;
  X=pinv(U)*temp;
  res=(-X(1,1)*(exp(K*Z)+exp(-K*Z))-(X(2,1)/(K*H1))*(1+K*Z)*exp(K*Z)+(X(3,1)/(K*H1))*(1-K*Z)*exp(-K*Z))*K*besselj(1,K*R);
endfunction

Je ne pense pas que le problème vienne de la fonction. Je l'ai testé avec des valeurs arbitraires et l’exécution de la fonction ne renvoie pas d'erreur.

Ensuite je fais appel à cette fonction de mon programme principale de la façon suivante :

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
 
nb_pts=10;
if strcmp(mode_exe,'yes')
    r=linspace(0,1E-3,nb_pts);
    z=linspace(0,7E-4,nb_pts);
else
    r=linspace(0,1E-3,nb_pts);
    z=linspace(0,7E-4,nb_pts);
end
 
%% Calculation of velocity field
% meshgrid
[V,W]=meshgrid(r,z);
% Velocity field urr
 
for i=1:length(V)
  for j=1:length(W)
    urr_1(i,j)=quadgk(@(k)BC_r1(k,W(j),V(i)),0,inf);
  end
end

J'ai aussi essayé de faire le calcul sans les deux boucles for en remplaçant "*" par ".*" et "/" par "./" où il fallait:

Code :

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

urr_1=quadgk(@(k)BC_r1(k,W,V),0,inf);

A l'exécution, octave me renvoie le message (c'est le même avec ou sans les boucles for) suivant :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
 
坤error: horizontal dimensions mismatch (150x3 vs 1x1)
error: called from:
error:   Z:\Donnees_these_DR\Interface_libre\Simulations_Calculs\thermocop_matlab\BC_r1.m at line 46,
 column 5
error:    at line -1, column -1
error:    at line -1, column -1
error: evaluating argument list element number 1
error:   C:\Program Files (x86)\Octave\Octave-3.8.1\share\octave\3.8.1\m\general\quadgk.m at line 425
, column 5
error:   C:\Program Files (x86)\Octave\Octave-3.8.1\share\octave\3.8.1\m\general\quadgk.m at line 307
, column 20
error:   Z:\Donnees_these_DR\Interface_libre\Simulations_Calculs\thermocop_matlab\thermocap_flows_mod
el_DR.m at line 45, column 15

Du coup je suppose que le problème vient de l'appel de la fonction BC_r1 dans quadgk, mais je ne comprend pas pourquoi.

[Beltharion]

Je ne pense pas que le problème vienne de la fonction. Je l'ai testé avec des valeurs arbitraires et l’exécution de la fonction ne renvoie pas d'erreur.

Du coup je suppose que le problème vient de l'appel de la fonction BC_r1 dans quadgk, mais je ne comprend pas pourquoi.

En fait je me suis complètement fourvoyé... Le problème vient de la fonction BC_r1.

Dans cette fonction j'utilise la fonction load (L24) pour faire appel à une structure contenant les valeurs des constantes.


Comment puis-je faire pour charger un fichier à l'intérieur de la fonction BC_r1? Et surtout est-ce que c'est possible ?

EDIT :
1) Le problème du chargement est réglé. L'erreur venait fichier que j'appelle qui n'avait pas le bon format

2) J'ai continué de chercher et j'ai identifié un autre problème (en espérant que ce soit le dernier...)
Il vient des matrices M et N définies dans la fonction BC_r1

Je les remet pour illustrer mon propos.

Code :

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

1
2
3
4
5
6
7
8
9
10
11
 
M=[-2 -1./(K*H1) -1./(K*H1) 2 1./(K*H2) 1./(K*H2); ...
    (exp(-K*H1)-exp(K*H1)) -exp(-K*H1) -exp(K*H1) 0 0 0; ...
    0 0 0 (exp(K*H2)-exp(-K*H2)) exp(K*H2) exp(-K*H2); ...
    (-exp(-K*H1)-exp(K*H1)) -((1-K*H1)/(K*H1))*exp(-K*H1) -((1+K*H1)/(K*H1))*exp(K*H1) 0 0 0; ...
    0 0 0 (-exp(K*H2)-exp(-K*H2)) -((1+K*H2)/(K*H2))*exp(K*H2) -((1-K*H2)/(K*H2))*exp(-K*H2); ...
    0 -2*(eta1/H1) 2*(eta1/H1) 0 2*(eta2/H2) -2*(eta2/H2)];
 
A=-sigma_vs_T.*K*Tmax*((wth.^2)/2).*exp(-((wth.*K).^2)/4);
[Li,Co]=size(A) %=> Li=1 et Co=150  au lieu de 1 !!!!
N=[0;0;0;0;0;A];

En théorie M est une matrice 6x6 et N une matrice 6x1. Or le problème est que le nombre de colonnes change d'une ligne à l'autre, et ce pour les deux matrices. Je ne comprend pas pourquoi !!!!
M et N sont écrites en fonction de K qui est une variable de la fonction BC_r1. Chaque élément des deux matrices doit prendre une valeur scalaire et non vectorielle.

En fait je ne suis pas certain de la source de l'erreur. Je me demande si le problème ne vient pas finalement de calcul itératif que je fais dans avec les boucles for dans le programme principal.

EDIT 2:
Quand j'exécute le programme principal, et que j'affiche la taille de K, j'obtiens une matrice 1x150 ! Je ne sais pas d'où viennent ces 150 colonnes !!!!
En revanche si je force la valeur de K, en écrivant par exemple K=K(1) (c'est qui totalement arbitraire et certainement faux...), les matrices M et N ont les bonnes tailles.

Cela m’emmène au problème suivant (après exécution du programme principal)

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
 
??? Error using ==> quadgk>finalInputChecks at 488
Output of the function must be the same size as the input.
 
Error in ==> quadgk>evalFun at 365
            finalInputChecks(x,fx);
 
Error in ==> quadgk>f2 at 395
        [y,too_close] = evalFun(t2t);
 
Error in ==> quadgk>vadapt at 276
            [fx,too_close] = f(x);
 
Error in ==> quadgk at 223
    [q,errbnd] = vadapt(@f2,interval);
 
Error in ==> thermocap_flows_model_DR at 35
    urr_1(i,j)=quadgk(@(k) BC_r1(k,z(j),r(i)),1E-12,inf);

D'après ce que j'ai pu lire, il existe une fonction arrayfun qui pourrait régler ce problème, mais je ne comprend pas son utilisation.

[Invité]

Avec arrayfun, on en reviendrait au même soucis des fonctions anonymes. (Son utilisation est basée sur le même principe que cellfun [cf FAQ] hormis qu'elle s'applique à chaque élément d'un tableau classique.)
On peut faire simple ici avec une boucle parcourant les éléments de K :

Code :

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

1
2
3
4
5
6
res = zeros(size(K));
for i = 1:numel(K)
    M = [ . . . K(i) . . .];
    % ...
    res(i) = . . .
end

[Beltharion]

Grâce à tes suggestions, j'ai réussi à faire fonctionner mon programme correctement.


Mais finalement, mon plus gros problème venait du dimensionnement des variables. En effet, en voulant faire un calcul en absolu, l'écart entre important entre les valeurs numériques des constantes faisait diverger le calcul.

Donc en adimensionnant le problème physique, le calcul converge. En plus dans mon cas, le fait d'avoir adimensionné mes données m'a libéré du caractère impropre de la transformée de Hankel.


Le problème est donc résolu !

Encore un grand merci à Winjerome