Débutant MatLab Résolu/Fermé

Question

Bonjour,
ma question risque de faire doucement sourire les plus aguerris d'entre vous...
Je cherche à créer une matrice n*n dont chacun des coefficients est fonction d'une variable. Je sais que MATLAB ne fonctionne pas en symbolique alors voici comment je procède pour l'instant:
"for w=0:0.01:100;
for i=1:n
    for j=1:n
        dav(i,j)=(0.03*u10^2./(pi.*wind(i,1).*wind(j,1).*w)).*(((1200.*w/u10)^2).*(1+(1200.*w/5.9).^2).^(-4/3)).*exp(-(w./u10).*C.*abs(z(i,1)-z(j,1)));
    end
end
end"
Le poblème est que je récupère l'ensemble des résultats les uns après les autres et que je souhaiterai au pire obtenir un vecteur colonne pour chacun des coefficients de la matrice, au mieux une matrice qui tient compte de la dimension "w".
J'espère que j'ai été assez clair dans mes explications.
Merci pour votre aide

Sacabouffe · Accepted Answer

Salut
Déjà, je viens de me rendre compte qu'il y avait une petite coquille dans ce que j'avais écrit pour la première méthode.
J'avais pris n=10 pour les premiers tests et j'ai oublié de changer dans le code.
Du coup, c'est plutôt ça:
function dav = mat_dav(w)
global C;
global u10;
global z;
global wind;
global n;
dav=zeros(n,n); %% C'est mieux de donner à MatLab la taille de la matrice au début
                %%  des boucles, c'est un peu plus rapide
for ii=1:n
for jj=1:n
dav(ii,jj)=(0.03*u10^2./(pi.*wind(ii,1).*wind(jj,1).*w)).*((­(1200.*w/u10)^2).*(1+(1200.*w/5.9).^2).^(-4/3)).*exp(-(w./u1­0).*C.*abs(z(ii,1)-z(jj,1)));
end
end Pour tracer une fonction sous MatLab, il faut de toute façon que tu lui donnes un intervalle et une discrétisation (donc le vecteur w). Du coup, si t'as utilisé les méthode du tableau à 3 dimensions, si n=10 et que tu veux tracer la courbes des valeurs du coefficient (2,3) par exemple, tu tapes:
plot(w,squeeze(dav(2,3,:))); Mais tu peux utiliser la méthode 1 et demander à MatLab de te tracer la courbe d'une fonction anonyme, ce qui n'est pas toujours une bonne chose puisque tu ne choisis pas la discrétisation.
Dans ce cas écris plutôt un truc du genre
function dav = mat_dav(w,indices)
global C;
global u10;
global z;
global wind;
global n;
if nargin==1
    dav=zeros(n);
for ii=1:10
for jj=1:10
dav(ii,jj)=(0.03*u10^2./(pi.*wind(ii,1).*wind(jj,1).*w)).*(((1200.*w/u10)^2).*(1+(1200.*w/5.9).^2).^(-4/3)).*exp(-(w./u10).*C.*abs(z(ii,1)-z(jj,1)));
end
end
else
    ii=indices(1,1);jj=indices(1,2);
    dav=(0.03*u10^2./(pi.*wind(ii,1).*wind(jj,1).*w)).*(((1200.*w/u10)^2).*(1+(1200.*w/5.9).^2).^(-4/3)).*exp(-(w./u10).*C.*abs(z(ii,1)-z(jj,1)));
end  Si tu rentres un nombres d'arguments égal à 1 (juste w), ça te donne la matrice dav pour un w donné, et si tu donnes en plus un indices (sous forme d'un vecteur 1x2), ça te donne le coefficient de la matrice pour les indices donnés.
Je te propose ça parce que je vois pas comment demander à MatLab de tracer seulement la courbe des valeurs d'un coefficient donné pour une fonction anonyme à valeurs matricielles :-(
Et aussi parce que je sais pas si dans tes applications tu as besoin d'avoir la matrice dav en entier.

Donc ensuite... si tu veux tracer les valeurs du coefficient (2,3) (donc n>=3) de ta matrice dav,
*soit tu fais
dav = @(w) mat_dav(w,[2,3]);
ezplot(dav) *soit tu choisis des valeurs xmin, xmax, ymin, ymax pour le graphique et tu tapes
dav = @(w) mat_dav(w,[2,3]);
fplot(dav,[xmin xmax ymin ymax]) Le @ sous MatLab, c'est pour définir les fonctions anonymes. Je pense que c'est ce que tu cherches quand tu parles de symbolique. Si tu veux plus d'infos, tu tapes help function_handle ou doc function_handle sous MatLab.

Voilà... tu me diras si ça marche...

Sacabouffe · Answer

Salut jubile_29
Je pense pas qu'il y ait à sourire, on a tous débuté en tout. En plus ta question est plutôt intéressante.
Juste un petit truc avant de commencer, évite d'utiliser i et j. Sous MatLab ce sont les i et j complexes. Dans ton programme ça pose pas de souci mais si un jour tu dois utiliser à la fois le i entier pour la boucle et le i complexe, tu risques d'avoir quelques soucis (et idem pour le j) ;-) 


1) Ce que tu peux faire, c'est définir une fonction en MatLab par exemple
function dav = mat_dav(w)
global C;
global u10;
global z;
global wind;
global n;
for ii=1:10
for jj=1:10
dav(ii,jj)=(0.03*u10^2./(pi.*wind(ii,1).*wind(jj,1).*w)).*(((1200.*w/u10)^2).*(1+(1200.*w/5.9).^2).^(-4/3)).*exp(-(w./u10).*C.*abs(z(ii,1)-z(jj,1)));
end
end global, c'est pour dire que les variables sont pas locales à la fonction, qu'elles sont définies ailleurs et que la fonction doit les utiliser. J'ai mis ça parce que je sais pas trop quelles sont ces variables. Si tu fais comme ça, il faudra à l'extérieur de la fonction que tu dises à MatLab que toutes ces variables sont globales et qu'ensuite tu leur attribues une valeur.

Si tu veux pas faire comme ça, tu peux tout simplement attribuer les valeurs qu'il faut à ces variables à l'intérieur de la fonction. Dans ce cas, vire toutes les lignes qui  sont là pour dire à la fonction qu'elles sont globales.

2) Il y a une autre solution, si tu veux pas de fonction mais directement les valeurs pour w=(0:0.01:100), tu peux stocker ça sous forme de tableau à trois entrées ou sous forme de cellules.

Tableau à 3 entrées:
w=(0:0.01:100);
for ii=1:n
for jj=1:n
for kk=1:size(w,2)
dav(ii,jj,kk)=(0.03*u10^2./(pi.*wind(ii,1).*wind(jj,1).*w(kk))).*(((1200.*w(kk)/u10)^2).*(1+(1200.*w(kk)/5.9).^2).^(-4/3)).*exp(-(w(kk)./u10).*C.*abs(z(ii,1)-z(jj,1)));
end
end
end Et là dav(:,:,7) par exemple, ça te donne la valeur de ta matrice pour w(7).

Cellules de matrices nxn:
w=(0:0.01:100);
for ii=1:n
for jj=1:n
for kk=1:size(w,2)
dav{kk}(ii,jj)=(0.03*u10^2./(pi.*wind(ii,1).*wind(jj,1).*w(kk))).*(((1200.*w(kk)/u10)^2).*(1+(1200.*w(kk)/5.9).^2).^(-4/3)).*exp(-(w(kk)./u10).*C.*abs(z(ii,1)-z(jj,1)));
end
end
end Et là dav{7} par exemple, ça te donne la valeur de ta matrice pour w(7).

Voilà, j'ai pas d'autres idées. J'espère que ça répond à peu près à ta question.

Par contre... les boucles en MatLab, c'est pas toujours génial, c'est lent... ici en particulier, 3 boucles emboîtées avec 10001 valeurs pour w, j'espère que ton n est pas trop grand, sinon ça doit patauger dans la semoule :-D
J'ai pas trop regardé comment on pouvait faire pour tes quelques lignes mais si tu pouvais vectorialiser, ce calcul, je pense que ce serait pas plus mal ;-)

Voili, voilou...
Bon après-midi

Sacabouffe · Answer

Salut

Ah OK... je m'étais pas trop intéressé à la signification de toutes tes variables dans ton code...
T'avais donc choisi w pour le &#969; de la pulsation et f est ta fréquence...

1) Si tu veux prendre plus de points dans [0;0.1] et moins dans [0.1;100], tu changes juste ton vecteur w
wlow=linspace(0,0.1,nwlow);
whigh=linspace(0.1,100,nwhigh);
w=[wlow whigh(2:end)]; En prenant soin de choisir les nombres de points nwlow et nwhigh. Ensuite tu lances le code pour calculer dav avec ce vecteur w et pour tracer le coefficient (2,3) de la matrice dav, tu fais comme précédemment
plot(w,squeeze(dav(2,3,:))); 2) Pour ton problème d'échelle logarithmique, tu te fixes une pulsation minimale wmin (pour exclure le 0), une pulsation maximale wmax (que t'avais prise égale à 100 au départ) et le nombre lw de plusations. Ensuite w=10.^(linspace(log10(wmin),log10(wmax),lw)); 3) Pour ce qui est de ton code qui patine un peu dans la semoule, comme je te disais, tu peux le vectorialiser. Pour les programmes avec de nombreuses boucles ça tourne généralement plus vite quand tu le fais. Par contre tu peux te heurter très vite à des problèmes de mémoire. Si tu veux avoir une idée de la manière dont s'y prendre, tu peux jeter un oeil ici:
http://www.commentcamarche.net/faq/sujet 11669 la vectorialisation sous matlab
Dans ton programme, je suppose que wind et z sont des vecteurs colonnes vu que seules la colonne 1 est utilisée. Du coup, j'écrirais ton programme vectorialisé de cette manière (ici j'ai pris l'échelle logarithmique pour w vu que je sais pas trop ce que tu vas faire)
w=10.^(linspace(log10(wmin),log10(wmax),lw));
wmat=permute(repmat(w,[n 1 n]),[3 1 2]);
windmat=repmat(wind*wind.',[1 1 sw]);
zmat=repmat(repmat(z,1,n)-repmat(z.',n,1),[1 1 sw]);
dav=(0.03*u10^2./(pi*windmat.*wmat)).*...
    (((1200*wmat/u10).^2).*(1+(1200*wmat/5.9).^2).^(-4/3)).*...
    exp(-(wmat/u10)*C.*abs(zmat)); Voilà... si t'as d'autres questions...

jubile · Answer

En fait, ce que je cherche à calculer c'est ce qui est présenté p.9 de l'article suivant:
https://web.mscsoftware.com/support/library/conf/auc99/p01099.pdf 
Jette un coup d'oeil si tu as un peu de temps je pense que ça vaut mieux que tous mes longs discours... Il s'agit en effet d'une intégration suivant les fréquences dont le résultat, si je comprend bien, est donné pour chacune des fréquences, au fur et à mesure de l'intégration. Je pense que c'est comme cela que ça fonctinne...

Sacabouffe · Answer

Ah OK... J'ai lu en diagonal...
C'est pas seulement la CRMS totale sur la plage de fréquence entière dont tu as besoin. Pour chaque fréquence f de la plage, il te faut la CRMS sur [0,f] si j'ai bien compris.
Du coup, à une constante additive près (et une constante multiplicative dont la valeur est le pas de fréquence), c'est plus un truc du genre
crms=sqrt(cumsum(PSD,3)); qu'il te faut.
Mais comme tout à l'heure, si t'as pris une échelle logarithmique pour les fréquences, faut être plus prudent.
C'est ce que tu cherches à faire ou pas?

jubile · Answer

Ben je ne sais pas si tu as bon parce que je n'ai pas encore testé l'opérateur cumsum mais en tout cas sur la théorie tu as compris l'essentiel. Je n'ai pas utilisé d'échelle logarithmique, mais plutôt:
"flow=linspace(0,0.1,1001);
 fhigh=linspace(0.1,100,1001);
 f=[flow fhigh(2:end)];"
qui st relativement efficace...
Je te tiens au courant des résultats que j'obtiens mais je pense que ça va merder parce que j'ai eu la bonne idée de multiplier ma matrice PSD (étape suivante de dav...) par une matrice complexe et sa conjuguée...inch allah, on verra...

najib_ingenieur · Answer

%racine carré de x
clc;
clear all;
x=input('donner un entier');
A=x;
e=1/2 * (x/A);
while (abs(e)>0.000001)
    A=A+e;
    e=1/2 * (x/A);
end
disp(['le racine carré de',num2str(x),'est',num2str(A)]);
end

Débutant MatLab

7 réponses

Discussions similaires