Discussion :

[domiq44]

Bonjour,

Je cherche à implémenter la fonction de Bessel "bessely()".

J'ai déjà implémenté la fonction "besselj()" comme ceci en m'inspirant des informations prises ici.

Voici le contenu du M-File "besselj.m":

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
 
function s=besselj(nu,z,nmax)
  for m=0:nmax-1
    a=(-1)^m;
    b=(z/2)^(2*m+nu);
    c=factorial(m);
    d=gamma(m+nu+1);
    e=a*b/c/d;
    J(m+1)=e;
  end
  s=sum(J);
end

Je l'ai testé en faisant ceci:

Code :

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

1
2
3
4
5
 
format long;
for z=0:.2:1
  besselj(1,z,100)
end

Et j'obtiens cela:

Code :

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

1
2
3
4
5
6
7
 
s = 0.0
s = 0.09950083263923601
s = 0.19602657795531875
s = 0.2867009880639158
s = 0.36884204609417
s = 0.44005058574493355

à comparer à ceci.

Est-ce que quelqu'un a déjà essayé d'implémenter la fonction "bessely()"?
Merci.

[domiq44]

Bonsoir,

J'ai trouvé une solution en utilisant une autre approche.

Je passe par les fonctions sphériques de Bessel jn(z) et yn(z), sachant que :

Code :

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

1
2
 
spherical bessel (n,z) = besselj(n+1/2,z)*sqrt(pi/(2*z))

Voici la solution pour les jn(z):

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
 
function y=jn(nu,z)
  if (nu < 0)
    y=-jn(-nu,z);
    return;
  end
 
  syms zz;
  a=(-zz)^nu;
  g=sin(zz)/zz;
  for m=1:nu
    g=diff(g,zz)/zz;
  end
  y=a*g;
  y=subst(z,zz,y);
end

Et celle pour les yn(z):

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
 
function y=yn(nu,z)
  if (nu < 0)
    y=-yn(-nu,z);
    return;
  end
 
  syms zz;
  a=-(-zz)^nu;
  g=cos(zz)/zz;
  for m=1:nu
    g=diff(g,zz)/zz;
  end
  y=a*g;
  y=subst(z,zz,y);
end

Si quelqu'un a d'autres solutions, ou des améliorations, je suis preneur.

Bonne soirée.

[Jerome Briot]

Dans quel but souhaites-tu réinventer la roue ?

Si tu as quelques compétences en Fortran, tu peux très bien utiliser la bibliothèque amos à la base des fonctions MATLAB.