appel d'une fonction

je reviens vers vous aujourd'hui avec un algorithme très complexe mais qui s’avère pertinent aussi. Avant de tout commencer voici mon algorithme:

Code:

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
m=10;
T=15;
U=15;
F=zeros(m+1,m+1);
 
for j=1:m
    for k=1:m
        F(j+1,k+1)=dblquad(@distribution,0,j*T/m,0,k*U/m);
    end
end
F1=F;
 
x=0:T/m:T;
y=0:U/m:U;
[xx,yy] = ndgrid(x,y);
knotsx = augknt([0:T/m:T],3)+T/(3*m);
knotsy = augknt([0:U/m:U],3)+U/(3*m);
 
zx=(2*x-T)/T;
 
syms z;
for i=1:m+1;
    w(i)=double(int(lagrange(z,i,zx),z,-1,1));
end
 
M=F;
 
 
 
Fn(m+1,m+1)=1;
 
iter=0;
while Fn(m+1,m+1)>0.0001
 
    iter=iter+1;
    Fn=zeros(m+1,m+1);
    for j=1:m+1
        for k=1:m+1
 
            xv2=((j-1)*T/m)-x;
            yv2=((k-1)*U/m)-y;
            [xxv2,yyv2] = ndgrid(xv2,yv2);
            d=distribution(xxv2,yyv2);
 
            sp = spap2(knotsy,4,y,F);
            coefsy = fnbrk(sp,'c');
            sp2 = spap2(knotsx,4,x,coefsy.');
            coefs = fnbrk(sp2,'c').';
            F=spcol(knotsx,4,x)*coefs*spcol(knotsy,4,y).';
 
 
            fonction=F.*d;
 
            for p=1:j
                for q=1:k
                    Fn(j,k)=Fn(j,k)+ ((j-1)*T/m) *((k-1)*U/m) *(1/4)*(mean(w)^2)*fonction(p,q);
                    %Fn(j,k)=Fn(j,k)+ ((j-1)*T/m) *((k-1)*U/m) *(1/4)*w(p)*w(q)*fonction(p,q);
                end
            end
        end
    end
    Fn=abs(Fn);
    F=Fn;
    M=M+F;
end
mesh(xx,yy,M);

avec le sous programme distribution contient cette fonction

Code:

1
2
3
4
5
6
7
m1=2;
    s1=0.5;
    m2=4;
    s2=1;
    rho=0.3;
 
    f=(1/(2*pi*s1*s2*sqrt(1-rho^2))) * exp( (-1/(2*(1-rho^2))) * ( ((x-m1).^2)/(s1^2) - 2*rho*(x-m1).*(y-m2)/(s1*s2) + ((y-m2).^2)/(s2^2)));

ce que je demande comme aide c'est de comment crée un sous programme à la place de 'distribution' et qui contient cette formulation

Code:

1
2
3
f=(1/(2*pi*s1*s2*sqrt(1-rho^2))) * exp( (-1/(2*(1-rho^2))) * ( ((x-m1).^2)/(s1^2) - 2*rho*(x-m1).*(y-m2)/(s1*s2) + ((y-m2).^2)/(s2^2)));
fi=inline(f);
ft=inline((1/((a^(o-1))^2))*fi(x/(a^(o-1)),y/(a^(o-1))));

avec la variable 'o' qui va incrémenter d'un cran à chaque appel, cela se fait automatiquement
merci pour vos aides