Discussion :

[follichon]

Bonjour à tous,

Voici un code qui j'ai écrit pour le calcul d'intégrale "à la Riemann". Ce n'est peut-être pas la meilleure façon de faire une intégrale double mais si quelqu'un a une suggestion, ça m'intéresse.

Merci

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
60
61
62
63
64
65
66
67
68
69
clear all
close all
t1=clock
pasxy=0.025;
[y,x]=meshgrid(-5:pasxy:5);
 
 
wo=4.25;
l=2.2;
k=2*pi/l;
zR=pi*wo^2/l;
xo=0;
zo=0;
a=0;
 
z=zo;
R=z*(1+(zR/z)^2);
w=wo*sqrt(1+(z/zR)^2);
P=atan(z/zR);
 
 
 
 
vx=[];
inf=-5;
sup=5;
pas=0.1;
 
inte=[];
 
for y=inf:pas:sup
 
 
    for x=inf:pas:sup
 
        X=(x+xo).*cos(a);
        Z=(x+xo).*sin(a);
 
        if isnan(R)
            Expo1=1;
            Expo2=1;
 
        else
            Expo1=exp(-i*k.*(X.^2+y.^2)./(2.*R));
            Expo2=exp(-i*k.*(x.^2+y.^2)./(2.*R));
 
        end
        clock
 
 
        A=sqrt(2/pi)*(1/w)*exp(-(X.^2+y.^2)/w.^2).*Expo1.*exp(-i*k*z+i.*P);
 
        Aprime=conj(A);
 
        B=sqrt(2/pi).*(1./w).*exp(-(x.^2+y.^2)./w.^2).*Expo2.*exp(-i*k*z+i.*P);
 
        C=Aprime.*B;
 
 
        vx=[vx x];
 
        inte1=sum(sum(abs(C)));
 
        inte=[inte inte1];
 
    end
end
 
t2=clock;

[Jerome Briot]

Quel est le problème ?
Quelle est la question exactement ?

[ol9245]

Tu peux supprimer toutes tes boucles. ca ira plus vite. Il te faut pour cela lancer les calculs sur le tableau entier et pas sur ses éléments un à un.