MonteCarlo Integrasjon

Det er god trening å prate matematikk. Her er det fritt fram for alle. Obs: Ikke spør om hjelp til oppgaver i dette underforumet.

Moderatorer: Vektormannen, espen180, Aleks855, Solar Plexsus, Gustav, Nebuchadnezzar, Janhaa

Svar
Nebuchadnezzar
Fibonacci
Fibonacci
Innlegg: 5648
Registrert: 24/05-2009 14:16
Sted: NTNU

Leste et artig bloggerinnlegg fra Fredrik angående montecarlo integrasjon.
Tenkte jeg kunne prøve å programmere dette selv.

Fikk det nesten til, men jeg får ikke plottet punktene mine. Kan noen si hvor feilen min ligger?

Og ja dette er Matlab =) Takker på forhånd for alle tips og hjelp

Kode: Velg alt

function MonteCarl( a )

% Genererer a antall x og y verdier mellom -1 og -1. 
x = 2*rand(a,1)-1;
y = 2*rand(a,1)-1;

% Setter sammen x og y koordinatene til punkter. 
h = [x y];
z = 0;
w = 0;

% Sorterer punktene som ligger utenfor enhetssirkelen og innenfor (r<1)
% Altså x^2 + y^2 < 1.
for i=1:a
    r = h(i,1)^2+h(i,2)^2;
    if r<1
        z = z + 1;
        p(z,1:2) = [ h(i,1) h(i,2) ];
    else
        w = w + 1;
        q(w,1:2) = [ h(i,1) h(i,2) ];
    end
end
% A er arealet av sirkelen. 
A = (z/a)*4; 
% D er avviket fra den egentlige verdien
D = abs(pi-A);

fprintf('Antall punkter innenfor sirkelen %d \n',z)
fprintf('Antall punkter utenfor  sirkelen %d \n',w)
fprintf('Sirkelens areal er ca %d med et avvik på %f  \n \n',A,D)

% Tegner en sirkel med radius 1. 
 [x,y,z] = cylinder(1,10^3);
 plot(x(1,:),y(1,:))
 hold on
 axis equal
 
 % Sjekker om vi har mer enn null punkter innenfor sirkelen
 % Om dette stemmer plotter vi disse punktene i rød.
 if z>0
 plot(p(:,1),p(:,2),'r+')
 end
 
 % Sjekker om vi har mer enn null punkter utenfor sirkelen 
 % dersom vi har det, plottes disse fargene i blå.
 if w>0
 plot(q(:,1),q(:,2),'b+')
 end
 hold off

end
"Å vite hva man ikke vet er og en slags allvitenhet" - Piet Hein
https://s.ntnu.no/Integralkokeboken
Lektor - Matematikk, Fysikk og Informatikk
drgz
Fermat
Fermat
Innlegg: 757
Registrert: 24/12-2008 23:22

Hva mener du med at du ikke får plottet punktene dine? Uten å endre koden får jeg opp sirkelen med haugevis av punkter som ligger på og utenfor sirkelen.

Ellers ser jeg at du overskriver variabelen z når du kaller
[x,y,z] = cylinder(1,10^3);

så derfor du sikkert ikke får noen røde punkter med koden din.

Ellers kan du sikkert skrive om til

Kode: Velg alt

r = h(:,1).^2+h(:,2).^2;
i1 = r<1;
i2 = ~i1;
zz = length(sum(i1));
w = length(sum(i2));
p = [h(i1,1) h(i1,2)];
q = [h(i2,1) h(i2,2)];
Sist redigert av drgz den 01/04-2012 16:13, redigert 1 gang totalt.
Nebuchadnezzar
Fibonacci
Fibonacci
Innlegg: 5648
Registrert: 24/05-2009 14:16
Sted: NTNU

Fant feilen min, min matlab likte vist ikke at jeg brukte z flere ganger med ulike betydninger. Om jeg forandrer z til s så fungerer alt

=)

EDIT:

Den kodesnippen du la til fungerte ikke, men en liten modifikasjon fikset problemet.

Kode: Velg alt

r = h(:,1).^2+h(:,2).^2;
i1 = r<=1;
i2 = r>1;
s = sum(i1);
w = sum(i2);
p = [h(i1,1) h(i1,2)];
q = [h(i2,1) h(i2,2)];
Må nok få bort min leie uvane med å elske løkker. Koden ble betydelig raskere å kjøre nå. Takker =D

En ganske stilig måte å beregne arealet av en sirkel på må jeg si.
"Å vite hva man ikke vet er og en slags allvitenhet" - Piet Hein
https://s.ntnu.no/Integralkokeboken
Lektor - Matematikk, Fysikk og Informatikk
drgz
Fermat
Fermat
Innlegg: 757
Registrert: 24/12-2008 23:22

Ja, ble noe tull her. Hadde først en kode med length() og ikke sum(), men så tenkte jeg at det ble bedre å bruke sum(), men glemte å ta vekk length() da jeg limte inn her, men du skjønte i alle fall hva jeg ville fram til. :)

Ellers ingenting i veien med løkker generelt (feks i andre språk), men i akkurat i matlab så er det en fordel å utnytte vektorer.
Nebuchadnezzar
Fibonacci
Fibonacci
Innlegg: 5648
Registrert: 24/05-2009 14:16
Sted: NTNU

Liten ting ting, matlab foretrekker

Kode: Velg alt

 [a1 a2] = size(p);
 plot(p(1:a1,1),p(1:a1,2),'r+')
over

Kode: Velg alt

plot(p(:,1),p(:,2),'r+')
Er det noen nevneverdig forskjell på disse to ?
"Å vite hva man ikke vet er og en slags allvitenhet" - Piet Hein
https://s.ntnu.no/Integralkokeboken
Lektor - Matematikk, Fysikk og Informatikk
FredrikM
Poincare
Poincare
Innlegg: 1367
Registrert: 28/08-2007 20:39
Sted: Oslo
Kontakt:

Artig at folk faktisk leser bloggpostene jeg skriver :D
En ganske stilig måte å beregne arealet av en sirkel på må jeg si.
Ikke veldig effektivt for sirkler, men ofte er det en veldig smart metode for kompliserte integrasjonsområder (en sirkel er vel noe av det minst kompliserte du kommer over).

Les f.eks de første avsnittene her: http://en.wikipedia.org/wiki/Monte_Carlo_methods

-

Kan (nesten) ingenting om Matlab, så kan ikke kommentere koden din. Selv bruker jeg som regel Sage (www.sagemath.org) når jeg skal visualisere. Gratis og open source (og at det er open source har jeg faktisk fått bruk for ;).
Cube - mathematical prethoughts | @MatematikkFakta
Med forbehold om tullete feil. (både her og ellers)
drgz
Fermat
Fermat
Innlegg: 757
Registrert: 24/12-2008 23:22

Nebuchadnezzar skrev:Liten ting ting, matlab foretrekker

Kode: Velg alt

 [a1 a2] = size(p);
 plot(p(1:a1,1),p(1:a1,2),'r+')
over

Kode: Velg alt

plot(p(:,1),p(:,2),'r+')
Er det noen nevneverdig forskjell på disse to ?
Hva legger du i "foretrekker"?
Nebuchadnezzar
Fibonacci
Fibonacci
Innlegg: 5648
Registrert: 24/05-2009 14:16
Sted: NTNU

Den gir meg en gul boks som sier at p forandrer lengde under hver iterasjon, og sier det er mer effektivt om jeg benytter meg av vektorer med fast lengde.

Så lurte jeg bare på om kjøretiden på disse to var den samme, eller om det er raskere å først beregne lengden.

Og ja, Fredrik jeg vet at det ikke er en spesielt effektiv metode å beregne arealet av en sirkel på, men det var en god oppfriskning i enkel programering. Men dog stilig =)
"Å vite hva man ikke vet er og en slags allvitenhet" - Piet Hein
https://s.ntnu.no/Integralkokeboken
Lektor - Matematikk, Fysikk og Informatikk
drgz
Fermat
Fermat
Innlegg: 757
Registrert: 24/12-2008 23:22

Nebuchadnezzar skrev:Den gir meg en gul boks som sier at p forandrer lengde under hver iterasjon, og sier det er mer effektivt om jeg benytter meg av vektorer med fast lengde.

Så lurte jeg bare på om kjøretiden på disse to var den samme, eller om det er raskere å først beregne lengden.
Hvis du bruker samme kode med for-løkke som i første posten så er det jo ikke rart. Du bruker p(z,1:2) = etellerannet uten å ha preallokert minne til matrisen p, så da må Matlab kopiere minne fram og tilbake for hver iterasjon, som ikke er så veldig praktisk. Hvis du har en annen kode nå er det greit å se på den først.
Nebuchadnezzar
Fibonacci
Fibonacci
Innlegg: 5648
Registrert: 24/05-2009 14:16
Sted: NTNU

Nå bruker jeg selvsagt ikke for løkker i koden min lengre (takk igjen), men jeg lurte på om det var noe forskjell i å først finne lengden av en vektor eller bare bruke :

Koden under brukte jeg for å teste dette.

Kode: Velg alt

x = 2*rand(10^6,1)-1;
y = 2*rand(10^6,1)-1;

h = [x y];

r = x.^2 + y.^2;
i1 = r<=1;
i2 = r>1;
s = sum(i1);
w = sum(i2);

P = [ h(i1,1) h(i1,2) ];
Q = [ h(i2,1) h(i2,2) ];

tic
[ n m ] = size(P);
P(1:n,1:2);
toc

tic
P(:,1:2);
toc

fprintf('\n')
Resultatet ble


Elapsed time is 0.015699 seconds.
Elapsed time is 0.009433 seconds.

Elapsed time is 0.014911 seconds.
Elapsed time is 0.008405 seconds.

Elapsed time is 0.015404 seconds.
Elapsed time is 0.008482 seconds.

Elapsed time is 0.016060 seconds.
Elapsed time is 0.008889 seconds.

Elapsed time is 0.016166 seconds.
Elapsed time is 0.008936 seconds.

Så det virker som det er litt raskere å benytte seg av :

I min kode valgte jeg heller å bruke

Kode: Velg alt

r = x.^2 + y.^2;
i1 = r<=1;
i2 = r>1;
s = sum(i1);
w = sum(i2);
og

Kode: Velg alt

 plot(h(i1,1),h(i1,2),'r*')
Kommer nok litt flere spørsmål fremmover angående matlab, fikk en ganske stor innlevering som må inn. Og virker som du har grei peiling. Der går det mer på nummerisk integrasjon, og plotting av figurer. Og mindre på løkker, vektorer og slikt. Men tenkte å fikle litt med det selv =)
"Å vite hva man ikke vet er og en slags allvitenhet" - Piet Hein
https://s.ntnu.no/Integralkokeboken
Lektor - Matematikk, Fysikk og Informatikk
drgz
Fermat
Fermat
Innlegg: 757
Registrert: 24/12-2008 23:22

Altså, greia var at du sier at Matlab klagde over at variabelen p endret størrelse for hver iterasjon; det betyr at du hadde en kode med løkker - man får ingen meldinger om sånt hvis man ikke bruker løkker. Så det var derfor jeg lurte på hva slags kode du egentlig brukte.

Ellers vil det som regel ta kortere tid hvis man bruker færrest mulig operasjoner (noe avhengig av hvilke operasjoner selvsagt).
Svar