CommentCaMarche
Recherche
Posez votre question Signaler

Débutant MatLab [Résolu]

jubile_29 4Messages postés mercredi 6 août 2008Date d'inscription 8 août 2008Dernière intervention - Dernière réponse le 22 oct. 2010 à 15:37
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
Lire la suite 
Réponse
+2
moins plus
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...
jubile_29 4Messages postés mercredi 6 août 2008Date d'inscription 8 août 2008Dernière intervention - 7 août 2008 à 14:21
Merci pour le "squeeze". C'est ce qu'il me manquait pour obtenir les courbes correspondant à chacun des coefficients de ma matrice en fonction de w. J'ai utilisé cette méthode qui est plus simple à mettre en application et à en comprendre les rouages... Je vais enfin pouvoir continuer ma programmation car, comme je te l'ai déjà dit, ce n'est malheureusement qu'une étape...
En tout cas, tu m'as enlevé une bonne épine du pieds
Merci encore mais il n'est pas impossible que j'ai encore besoin de tes lumières par la suite (sans vouloir abuser de ton temps).
Bonne aprèm
Répondre
Sacabouffe 9603Messages postés dimanche 19 août 2007Date d'inscription 29 mai 2009Dernière intervention - 7 août 2008 à 14:24
Pas de souci ;-)
Si c'est une question un peu différente, pense à ouvrir un nouveau topic, ce sera pas plus mal. Et puis, pour une autre question, il est pas impossible que j'ai pas la réponse, en ouvrant un nouveau sujet t'auras plus de chances d'avoir l'aide d'une personne qui pourra te conseiller correctement.
Bonne journée !
Répondre
jubile_29 4Messages postés mercredi 6 août 2008Date d'inscription 8 août 2008Dernière intervention - 8 août 2008 à 14:43
Salut,
j'ai revérifié mon code MATLAB et je m'étais laissé berner par la véracité des résultats que j'obtenais...
Le plot(f,squeeze(dav(2,3,:))) ne fonctionne pas car les length de f et dav ne sont pas identiques. En revanche, lorsque j'utilise cette matrice à trois entrées dav pour définir
PSDin(ii,jj,kk)=4*v(ii,1).*v(jj,1).*dav(ii,jj,kk) où v(ii,1) est un vecteur 5*1
alors seulement dans ce cas PSDin est de la même longueur que f et je peux donc tracer mes courbes. C'est à n'y rien comprendre!!! En tout cas je suis perdu...et je me demande maintenant si mes calculs sont corrects d'autant plus que j'ai besoin d'utiliser le même cheminement pour construire une autre matrice à trois entrées mais que lorsque je souhaite tracer cette nouvelle matrice, je me retrouve coincé par le même problème qu'évoqué ci-dessus, à savoir un problème de longueur...
Merci pour votre aide
Répondre
Sacabouffe 9603Messages postés dimanche 19 août 2007Date d'inscription 29 mai 2009Dernière intervention - 8 août 2008 à 18:35
Salut
Il y a pas de f dans tout ce que m'as donné :-(
C'est quoi f?
Répondre
jubile- 11 août 2008 à 10:06
Ah oui, pardon...
En fait le "f" c'est le "w" de départ. J'ai tout simplement opéré nu changement de variable dans mes équations en effectuant "w=2pif". Ca doit te rappeler des choses, non? Sinon, j'ai réussi à trouver mon problème, qui s'est avéré être un réel problème de débutant, à savoir que j'oubliais de réinitialiser mon calcul à chaque fois que je le relançais etça, MATLAB ne l'aimait pas.
J'ai un autre petit soucis concernant la représentation de mes résultats. Comme tu l'avais prévu, ça patauge pas mal dans la semoule compte tenu du nombre de boucle "for" et même si j'obtient des résultats, ceux-ci ne sont pas complètement exploitable par ce que mon pas en "w/f" est trop grand pour les "w/f" très faibles (entre 0 et 0.1). Je souhaite prendre lpus de points dans cet intervalle. Ma première idée était de déterminer deux intervalles: le premier avec un pas de temps très petit entre 0 et 0.1 et le second avec un pas de temps plus grand entre 0.1 et 100, puis de tracer sur le même graphe mes deux courbes-résultats dav1 et dav2 correspondant aux deux intervalles. Malheureusement cela ne fonctionne pas (peut-être est ce du à mon encodage, d'ailleurs!)
Connaitrais-tu une solution simple afin de définir un tel pas de temps (du genre logarithmique...), c'est-à-dire très petit pour les petites valeurs de "w/f" et plus grands pour les valeurs supérieures...
Merci
Répondre
Ajouter un commentaire
Réponse
+0
moins plus
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
jubile_29 4Messages postés mercredi 6 août 2008Date d'inscription 8 août 2008Dernière intervention - 7 août 2008 à 11:53
Merci pour ton aide Sacabouffe,
je récupère bien la matrice qu'il faut avec la deuxième méthode que tu proposes. Cependant, je suis loin d'être sorti d'affaire puisque l'aboutissement de mon application physique est la construction d'un graphe fonction de w de chacun des coefficients de ma matrice généralisée en w (cette matrice dav n'est d'ailleurs qu'une partie de ma matrice finale en w que je souhaite calculer sous matlab). C'est pourquoi je suis si déçu qe MATLAB ne comprenne pas le symbolique et que je pensais que ta première idée(définition d'une fonction en w) pouvait convenir davantage à mon problème. Malheureusement, je n'arrive pas à la mettre en application car lorsque j'utilise directement la programmation que tu as décrite, en entrant directement les valeurs de mes constantes, MATLAB me demande de définir w or je pensais qu'il n'avait pas besoin de cela pour fonctionner...
Vois tu un moyen de récupérer chacun des coefficient de la matrice dav et de tracer un graphe de ce coefficent en fonction de w?
Merci encore
Répondre
Ajouter un commentaire
Réponse
+0
moins plus
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 ω 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- 11 août 2008 à 15:59
Wahou!!!ok meci, c'est nettement plus rapide lorsque j'utilise le linespacepour les points dont j'ai besoin. Etant donné que mes calculs ne sont pas très compliqués et longs de manière général (quelques dizaines de secondes en utilisant le bon échantillonnage)je n'ai pas particulièrement besoin de vectorialiser mon code...Du coup, je me sens un peu gêné que tu ais fourni tant d'effort d'autant plus que j'ai encore une question concernant ce code pour que j'obtienne, normalement, l'ensemble des résultats nécessaires. Je souhaite calculer la racine carrée de la somme en "w ou f" des coefficients de ma matrice. Et ensuite, plotter ces résultats en fonction de w ou f. Comment dois-je m'y prendre, en sachant que la fonction quad ne marche pas, que j'ai également essayé trapz et cumtrapz?...
Merci encore
Répondre
Sacabouffe 9603Messages postés dimanche 19 août 2007Date d'inscription 29 mai 2009Dernière intervention - 11 août 2008 à 18:00
De rien ;-)
Pour la suite, tu fais
rsdav=squeeze(sqrt(sum(sum(dav))));
plot(w,rsdav);
Ça devrait aller...
Répondre
jubile- 12 août 2008 à 10:26
Salut,
"rsdav=squeeze(sqrt(sum(sum(dav))));
plot(w,rsdav);" ne me permet pas, me semble t-il, d'obtenir les résultats pour chacun des coefficients mais sur l'ensemble de ces coefficients, non? Je ne sais pas dans quel branche tu travailles, ni même si tu as des notions de dynamique aléatoire ou traitement du signal, mais, en fait, mon dav est un élément de la PSD (Power of Spectral Density) de déplacement d'un système mécanique que j'ai déjà défini à l'aide de tes conseils. Chacun des coefficients de la diagonale de mes matrices en f obtenu représente donc un point de la PSD de déplacement d'un point précis du système étudié. Le fait de tracer en fonction de f (sur une certaine gamme) me permet donc d'obtenir la PSD de déplacement complète sur cette bande fréquentielle en un point de mon système. Le but recherché par la dernière question posée est de tracer la CRMS (Cumulativ Root Mean Square) de déplacement en fonction de f afin de déterminer l'évolution du déplacement-cette fois-ci-d'un point de mon système (donc d'un des coefficients de la diagonale) en fonction de la fréquence d'excitation. Je sais que MATLAB connaît les fonctions PSD et CRMS mais, dans ma logique et afin de bien comprendre le fonctionnement de mon système, je préfère rentrer moi-même un code équivalent qui me donne les résultats attendus. Quoiqu'il en soit, lorsque je fais:
"RMS=squeeze(sqrt(sum(sum(dav(5,5,:))))) %j'utilise ici dav afin que tu puisses comprendre le pb de manière général parce qu'en réalité, cette matrice-fonction n'est plus qu'un problème résolu%
loglog(f,RMS)"
de 1, le résultat obtenu pour RMS est unique alors que je m'attends à obtenir un vecteur colonne où le nombre d'éléments est égal à mon échantillonnage (cad nb de f)
et comme il trouve un résultat indéterminé, il ne peut pas me tracer la fonction (mais c'est un autre problème je pense)...
J'espère que tu as compris ce que je veux dire
Merci
Répondre
Sacabouffe 9603Messages postés dimanche 19 août 2007Date d'inscription 29 mai 2009Dernière intervention - 12 août 2008 à 13:44
Salut jubile
J'avais mal compris ta ta demande.

Je souhaite calculer la racine carrée de la somme en "w ou f" des coefficients de ma matrice.

J'avais compris que tu devais calculer pour chaque fréquence (mon interprétation de "en f"), la racine carrée de la somme des coefficients de la matrice.

Du coup, il y a un truc que j'ai un peu de mal à saisir. Soit tu fais une CRMS sur les déplacements et dans ce cas, il y a plus de dépendance en fonction des coefficients, soit tu fais une CRMS sur les fréquences (donc sur le temps d'après le théorème de Parseval) et dans ce cas, il y a plus de dépendance en fonction des fréquences.
Je comprends pas comment tu veux avoir à la fois une dépendance sur les coefficients de la matrice et une dépendance sur les fréquences.

En tout cas, si tu veux une CRMS sur les fréquences, c'est
crms=sqrt(sum(PSD,3));
Et là, t'as la CRMS pour tous les coefficients de ta matrice.

Si t'as pris une subdivision régulière pour les fréquences, c'est en fait la CRMS à un facteur près (la pas de fréquence).
Si t'as pris une subdivision logatithmique, c'est complètement faux. Faut adapter un peu pour l'intégation numérique.

Par contre, si c'est une CRMS sur les déplacements, je pense que c'est plutôt la racine carrée de la somme des coefficients au carré de la matrice, non?
Répondre
Ajouter un commentaire
Réponse
+0
moins plus
En fait, ce que je cherche à calculer c'est ce qui est présenté p.9 de l'article suivant:
http://www.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...
Ajouter un commentaire
Réponse
+0
moins plus
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?
Ajouter un commentaire
Réponse
+0
moins plus
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...
Sacabouffe 9603Messages postés dimanche 19 août 2007Date d'inscription 29 mai 2009Dernière intervention - 12 août 2008 à 17:33
Si t'as pris ça pour la discrétisation, il faut adapter aussi.
crms=sqrt(cumsum(cat(3,0.0001*PSD(:,:,1:1001),0.0999*PSD(:,:,1002:end)),3));
Ça devrait marcher...
Répondre
Ajouter un commentaire
Réponse
+0
moins plus
%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
Ajouter un commentaire
Ce document intitulé «  Débutant MatLab  » issu de CommentCaMarche (www.commentcamarche.net) est mis à disposition sous les termes de la licence Creative Commons. Vous pouvez copier, modifier des copies de cette page, dans les conditions fixées par la licence, tant que cette note apparaît clairement.

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes.

Le fait d'être membre vous permet d'avoir des options supplémentaires.