Ubuntu-fr

bochra017

Svp Aide code Matlab

Bonjour,

j'ai un code Matlab pour 1 dimension, j'ai besoin de votre aide pour écrit même code pour 2 dimension.

J’ai besoin de montrer que l'erreur Monte-Carlo est fonction de 1/sqrt(N) . Mon prof me dit que : Pour obtenir les erreurs MC, pour chaque nombre de points N (10^2, 10^3, 10^4...10^7) vous refaites le calcul (30 fois par exemple) pour calculer sigma, et l'erreur sera sigma/sqrt(N). vous aurez l'erreur pour chaque N, vous tracer la courbe pour obtenir la courbe en bleu. Et tracer 1/sqrt(N) qui est la courbe verte. Il me donne un exemple d’intégrale pour 1D ( code+commentaires explicatif ) et il me demande de faire la même chose pour 2D, pour montré que l’erreur converge vers 1/sqrt(N) n'importent qu’il dimension. Mais je n’arrive pas affaire , j’ai besoin l'aide svvvp.

code 1D :

close all
clc
format long

N = [10^2 10^3 10^4 10^5 10^6 10^7]';
C = 0.25; %fmax
a = 0;
b = 1;

%Nombre des essais
MCTries = 7;

%On aura besoin de ces 2 tableaux pour garder enregister lerreur 
%correspondante à chaque essaye (pour cet exemp: pour chaque N on refera le calcul 7 fois.
%on prendera à la fin la moyenne des erreurs.
array_sigmas = zeros(MCTries,1);
array_errors = zeros(length(N),1);

for i=1:length(N) % boucle sur les N
    for j=1:MCTries % boucle sur les essais
        I = 0;
        f = zeros(N(i,1),1);
        for k = 1:N(i,1) % boucle sur les xi
            x = rand();
            fx = x*(1-x)*sin(200*x*(1-x))^2;
            I = I + fx;
            f(i,1) = fx;%stockage des fx
        end
        % I = (b-a)*I/N(i,1) % valeur de l'integrale
    array_sigmas(j,1) = sqrt(var(f)); % calcul de l'erreur (la racine de la variance)  
    end
    % on prends la moyenne des essais
    array_errors(i,1) = mean(array_sigmas);
end


%data to plot
 x = N(:,1);
y0 = 1./sqrt(N(:,1));
y1 = array_errors(:,1)

%plotting
figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
box(axes1,'on');
grid(axes1,'on');
plot(x,y0,'-r',x,y1,'-bo');
set(axes1,'XMinorTick','on','XScale','log','YMinorTick','on','YScale','log');
xlabel('N') % x-axis label
ylabel('errors') % y-axis label
legend('1/sqrt(N)','error in MC integration')



<

pour 2D : par exemple cette fonction f(x,y)=xye^((-x^2)y) ou bien une fonction de votre choix

Et merci.

grigouille

Re : Svp Aide code Matlab

matlab n'est pas dans les dépôts donc, je ne sais pas si tu auras beaucoup d'aide.
Après, ce type de problème peut très bien se faire dans un autre language (C, C++, Fortran, etc.)

J'espère que tu es confiant sur ton générateur de nombres aléatoires.
Teste le d'abord s'il est vraiment bon sur autant de valeurs.

---

Debian (xfce) 12
HP LaserJet M1132 MFP

bochra017

Re : Svp Aide code Matlab

Il faut que je résoudre ce problème avec Matlab, car c'est  prof qui ma dit ça. c’est une partie de mon  projet . dit moi ou je vais aller pour avoir une aide, Est-ce que tu connais quelqu'un qui peut m’aider ?

bochra017

Re : Svp Aide code Matlab

J'ai vraiment besoin de votre aide svp! ,tout ce que je veux c'est :même code pour 2D.

DonutMan75

Re : Svp Aide code Matlab

Hello,
c'est un peu difficile de comprendre où tu bloques, tu as simplement copié/collé l'énoncé du TD...
Bloques-tu sur la compréhension "mathématique" du calcul d'une intégrale 2D ?
Bloques-tu sur l'implémentation Matlab ?
Qu'as-tu déjà tenté ?

Une première piste serait par exemple de partir de la ligne

x = rand();

et de la remplacer par :

x = rand();
y = rand();

Donut

Dernière modification par DonutMan75 (Le 29/06/2017, à 15:38)

bochra017

Re : Svp Aide code Matlab

j'ai bloque  sur l'implémentation Matlab, j'ai pas le temp d'apprendre les cours de Matlab.                    tout ce que je veux  c'est  d'avoir le figure pour 2D, j'ai besoin de votre aide pour  écrire le code en 2D.

ce qui concerne :

%data to plot
 x = N(:,1);
y0 = 1./sqrt(N(:,1));
y1 = array_errors(:,1)

cette partie qui ma fait tromper, je ne sais pas comment la remplacer.

et merci d'avance

DonutMan75

Re : Svp Aide code Matlab

Hello,
alors dans ce cas tu vois algorithme à développer mais ne sait pas comment l'implémenter en Matlab ?
Il faudrait donc que tu nous donnes le pseudo-code et on se fera un plaisir de t'aider à convertir ça en Matlab.
Par contre on va pas partir de l'énoncé copié/collé de ton TD et te le résoudre à ta place...

De ce que je comprends du machin, la grosse modif doit se passer dans la boucle for sur les xi

for k = 1:N(i,1) % boucle sur les xi
[...]

Si tu remplaces les xi par un point aléatoire 2D (comme je l'indiquais tantôt), n'est-ce pas un premier pas vers la solution ?

Donut

bochra017

Re : Svp Aide code Matlab

ok, je vais te donner l'algorithme de méthode , ce qui concerne l'erreur j'ai justement le code pour 1D .
algorithme de methode

pour 2D: je vais essayer de remplacer

C = 0.25; %fmax

par :

C = 3; %fmax

et

for i=1:length(N) % boucle sur les N
    for j=1:MCTries % boucle sur les essais
        I = 0;
        f = zeros(N(i,1),1);
        for k = 1:N(i,1) % boucle sur les xi
            x = rand();
            fx = x*(1-x)*sin(200*x*(1-x))^2;
            I = I + fx;
            f(i,1) = fx;%stockage des fx
        end

par exemple, par :

for i=1:length(N) % boucle sur les N
    for j=1:MCTries % boucle sur les essais
        I = 0;
        f = zeros(N(i,1),1);
        for k = 1:N(i,1) % boucle sur les xi et yi
for k = 1:N(i,1)
            x = rand();
            y = rand();
           fx = x.*y.*exp(((-x).^2).*y);
            I = I + fx;
            f(i,1) = fx;%stockage des fx
        end

et .

%data to plot
 x = N(:,1);

par:

%data to plot
 x = N(:,1);
 y = N(:,1);

mais je n'ai pas compris N(:,1) !!!!!

Dernière modification par bochra017 (Le 29/06/2017, à 21:20)

DonutMan75

Re : Svp Aide code Matlab

Hello,
De ce que je comprends N est un vecteur contenant les différents nombres de point pour les simulations.
N est donc un vecteur 1D.
Partant de là, toutes les références à N(i,1) me paraissent suspectes (bien que fonctionnelles) et pourraient avantageusement être remplacée par N(i) directement. Ca augmente la lisibilité du code.

Après cette remarque, une fois que tu as implementé la modif dans la boucle for k = 1:N(i,1), je ne pense pas que tu aies quoi que ce soit à changer dans le reste du code.
Est-ce que ça génère une erreur à l'exécution ?

Donut

bochra017

Re : Svp Aide code Matlab

Mercii bcp monsieur pour votre aide ,

close all
clc
format long

N = [10^2 10^3 10^4 10^5 10^6 10^7]';
C = 3; %fmax
a = 0;
b = 1;

%Nombre des essais
MCTries = 7;

%On aura besoin de ces 2 tableaux pour garder enregister lerreur 
%correspondante à chaque essaye (pour cet exemp: pour chaque N on refera le calcul 7 fois.
%on prendera à la fin la moyenne des erreurs.
array_sigmas = zeros(MCTries,1);
array_errors = zeros(length(N),1);

for i=1:length(N) % boucle sur les N
    for j=1:MCTries % boucle sur les essais
        I = 0;
        f = zeros(N(i),1);
        for k = 1:N(i) % boucle sur les xi et yi
            x = rand();
            y = rand();
            fx = x.*y.*exp(((-x).^2).*y);
            I = I + fx;
            f(i,1) = fx;%stockage des fx
        end
        % I = (b-a)*I/N(i,1) % valeur de l'integrale
    array_sigmas(j,1) = sqrt(var(f)); % calcul de l'erreur (la racine de la variance)  
    end
    % on prends la moyenne des essais
    array_errors(i,1) = mean(array_sigmas);
end


%data to plot
 x = N(:,1);
 y = N(:,1);
y0 = 1./sqrt(N(:,1));
y1 = array_errors(:,1)

%plotting
figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
box(axes1,'on');
grid(axes1,'on');
plot(x,y0,'-r',x,y1,'-bo');
set(axes1,'XMinorTick','on','XScale','log','YMinorTick','on','YScale','log');
xlabel('N') % x-axis label
ylabel('errors') % y-axis label
legend('1/sqrt(N)','error in MC integration')

https://up.harajgulf.com/imagef-1498832 … 1-png.html

j'obtient des très bon résultats, mais ce qui concerne 

plot(x,y0,'-r',x,y1,'-bo');

!!! cette partie je ne sais pas comment la remplacer.

et merci d'avance

Modération : merci à l'avenir d'utiliser les balises code (explications ici).

Dernière modification par cqfd93 (Le 03/07/2017, à 15:24)

DonutMan75

Re : Svp Aide code Matlab

Hello, j'ai regardé un peu plus en profondeur ton code et j'ai quelques remarques :

- tout d'abord l'exemple initial 1D, c'est ton prof qui te l'a donné ??? Il me paraît très bancal....

- la variable C n'est pas du tout utilisée dans le script initial, à quoi sert-elle ?

- le code ci-dessous est commenté... C'est pourtant lui qui calcule la valeur de l'intégrale par tirage aléatoire

%I = (b-a)*I/N(i,1)

I/N(i, 1) --> ça c'est la valeur moyenne du signal
(b-a)*I/N(i, 1) --> ça c'est l'estimation de l'intégrale

- f(i) est un vecteur de taille length(N), f(i) contient la dernière valeur tirée de façon alétoire (i.e. dernière itération de la boucle sur les k) ==> Il y a sans doute une erreur ici...
- array_sigmas(j) est un vecteur de taille MCTries, il contient la racine de la variance du signal f, mais ça n'a aucun sens ???
Si je comprends bien ce que tu cherches à faire, il faudrait plutôt :
f = vecteur d'erreur entre estimation de l'intégrale avec N(i) tirages aléatoire d'un côté, et intégrale "vraie" de l'autre (par exemple estimé avec méthode des rectangles avec petit pas de temps OU solution analytique)
array_sigma = les différents f selon les MCTries (dont array_sigma de taille MCTries)
et finalement array_errors de taille length(N) puisqu'on prend la moyenne pour chaque N(i)

Edit : OU ALORS calculer la variance en fonction du MCTries. Avec un N assez élevé, on retrouve la même valeur estimée d'un essai à l'autre... Alors qu'avec un N trop faible, on fait apparaître des différences notables... Ca évite de devoir comparer à un signal de référence.

Bref, je pense que tu bloques plus sur un problème mathématique que sur un problème de Matlab.

Essaie de découper le problème en plusieurs étapes.
Fixe tout d'abord N (le nombre de tirage aléatoire) à une valeur donnée (par exemple 100) et ne fais qu'un seul essai (MCTries=1) : comment estimer par tirage aléatoire la valeur de l'intégrale d'une fonction donnée et comment estimer l'erreur associée à cette estimation ?

Dans un second temps, on peut envisager de faire varier N et MCTries pour faire des estimations statistiques sur la qualité de la méthode.

Pour terminer, pourquoi ne pas estimer la valeur de l'intégrale par méthode des rectangles ? C'est moins coûteux en terme de calcul et la disposition des points sur le domaine de définition est plus homogène. Il ne reste que le côté didactique du truc, histoire de programmer un petit tirage aléatoire....

Enfin bref, bon courage pour la suite !

Donut

Dernière modification par DonutMan75 (Le 03/07/2017, à 09:58)