Einführung in Mathematica XVII

Integrieren mit Zufallszahlen

© Jürgen Schnack, 2010/05/29

Problem

Integrale mit kompliziertem Integraden bzw. in hohen Dimensionen kann man oft nicht analytisch berechnen. Oft ist "Würfeln" die Methode der Wahl! 
Weil diese Methoden auf Zufallszahlen basieren, heißen sie 
Monte-Carlo-Methoden.

1. Beispiel

f[x_]=1.0+(E^(Sin[x]));

Pf=Plot[f[x],{x,0,10},
PlotRange->{0,4},
PlotStyle->{Red,Thick},
BaseStyle->{FontFamily->"Helvetica",FontSize->18}]

Integrate[f[x],{x,0,10}]

Natürlich ginge NIntegrate, aber das wollen wir jetzt nicht.

Idee (hit or miss)

BG=Graphics[{Blue,Opacity[0.1],
Rectangle[{0,0},{10,4}]}];

LP1=ListPlot[Table[{RandomReal[{0,10}],RandomReal[{0,4}]},
{n,1,(10^2)}],PlotRange->{{0,10},{0,4}},
AspectRatio->1,BaseStyle->{FontFamily->"Helvetica",FontSize->18}];

Show[Pf,BG,LP1]

Punkte würfeln und auszählen!

Monte-Carlo-Integration (Hit or Miss)

NGes=10000;
NGood=0;
AGes=4.0*10.0;
For[i=0,i<NGes,i++,
xr=RandomReal[{0.0,10.0}];
yr=RandomReal[{0.0,4.0}];
(* Print[xr,"  ",yr,"  ",f[xr]]; *)
If[yr≤f[xr],NGood=NGood+1]]
AGood=(NGood/NGes)*AGes

NIntegrate[f[x],{x,0,10}]

Zweidimensionales Beispiel I

Use Green's theorem to find the area of a parametrically defined two-dimensional region:

x[t_]:=Cos[t]+1/7Cos[7t+Pi/3]; 
y[t_]:=Sin[t]+1/7Sin[7t]; 
ZDP1=ParametricPlot[{x[t], y[t]}, {t, 0, 2π}]

NIntegrate[x[t]D[y[t],t],{t,0,2π}]

Zweidimensionales Beispiel II

ZDBG=Graphics[{Blue,Opacity[0.1],
Rectangle[{-1.5,-1.5},{1.5,1.5}]}];

ZDLP1=ListPlot[Table[{RandomReal[{-1.5,1.5}],
RandomReal[{-1.5,1.5}]},
{n,1,(10^2)}],PlotRange->{{-1.5,1.5},{-1.5,1.5}},
AspectRatio->1,BaseStyle->{FontFamily->"Helvetica",FontSize->18}];

Show[ZDP1,ZDBG,ZDLP1]

Aber wie stellt man jetzt fest, ob ein Punkt innerhalb oder außerhalb liegt?
Ideen??? Für Lösungen gibt es "Echte-Programmierer-Punkte".

2. Beispiel

f2[x_]=Tan[x]^2;

Pf2=Plot[f2[x],{x,0,1.3},
PlotRange->{0,10},
PlotStyle->{Red,Thick},
BaseStyle->{FontFamily->"Helvetica",FontSize->18}]

Monte-Carlo-Integration (Mittelwertmethode)

NGes=10000;
FSum=0.0;
For[i=0,i<NGes,i++,
xr=RandomReal[{0.0,1.3}];
FSum=FSum+f2[xr]];
FInt=FSum*1.3/NGes

Integrate[f2[x],{x,0,1.3}]

Monte-Carlo-Integration: Varianz

jmax=5;
IntTable=Table[{(10^j),0.0},{j,1,jmax}];
VarTable=Table[{(10^j),0.0},{j,1,jmax}];
For[j=1,j<jmax+1,j++,
NGes=(10^j);
FSum=0.0;
FSumSq=0.0;
For[i=0,i<NGes,i++,
xr=RandomReal[{0.0,1.3}];
FSum=FSum+f2[xr];
FSumSq=FSumSq+(((f2[xr]))^2)];
FInt=FSum*1.3/NGes;
IntTable[[j,2]]=FInt;
FVar=(FSumSq/NGes-(((FSum/NGes))^2))/NGes;
VarTable[[j,2]]=(Sqrt[FVar]);
]

IntLP=ListLogLogPlot[IntTable,
PlotStyle->{Hue[0.0],PointSize[0.015]},BaseStyle->{FontFamily->"Helvetica",FontSize->18}]

VarLP=ListLogLogPlot[VarTable,
PlotRange->{{5,2*(10^jmax)},{0.001,1.0}},
PlotStyle->{Hue[0.0],PointSize[0.015]},BaseStyle->{FontFamily->"Helvetica",FontSize->18}]

2.3 ist das aber immer noch nicht!!!

Monte-Carlo-Integration (importance sampling)

w[x_]=2*(x^6);

Pf2=Plot[{f2[x],w[x]},{x,0,1.3},
PlotRange->{0,15},
PlotStyle->{Red,Thick},
BaseStyle->{FontFamily->"Helvetica",FontSize->18}]

Pf2=Plot[f2[x]/w[x],{x,0.01,1.3},
PlotRange->{0,15},
PlotStyle->{Red,Thick},
BaseStyle->{FontFamily->"Helvetica",FontSize->18}]

yw[x_]=Integrate[w[t],{t,0,x}]

xw[y_]=(((7/2)y)^(1/7));

NGes=10000;
FSum=0.0;
For[i=0,i<NGes,i++,
yr=RandomReal[{0.0,yw[1.3]}];
FSum=FSum+f2[xw[yr]]/w[xw[yr]]
];
FInt=FSum*(yw[1.3]-yw[0.0])/NGes

Integrate[f2[x],{x,0,1.3}]

Hm, so überzeugend ist das aber auch nicht!!!

Monte-Carlo-Integration (importance sampling-Varianz)

jmax = 5;
ISIntTable = Table[{10^j, 0.0}, {j, 1, jmax}];
ISVarTable = Table[{10^j, 0.0}, {j, 1, jmax}];
For[j = 1, j < jmax + 1, j++,
 NGes = 10^j;
 FSum = 0.0;
 FSumSq = 0.0;
 For[i = 0, i < NGes, i++,
  yr=RandomReal[{0.0,yw[1.3]}];
FSum=FSum+f2[xw[yr]]/w[xw[yr]];
FSumSq = FSumSq + (f2[xw[yr]]/w[xw[yr]])^2];
FInt=FSum*(yw[1.3]-yw[0.0])/NGes;
ISIntTable[[j, 2]] = FInt;
 FVar = (FSumSq/NGes - (FSum/NGes)^2)/NGes;
 ISVarTable[[j, 2]] = Sqrt[FVar]]

ISVarTable

IntLP = ListLogLogPlot[ISIntTable,
  PlotStyle -> {Hue[0.0], PointSize[0.015]}, 
  BaseStyle -> {FontFamily -> "Helvetica", FontSize -> 18}]

VarLP = ListLogLogPlot[ISVarTable,
  PlotRange -> {{5, 2*10^jmax}, {0.001, 1.0}},
  PlotStyle -> {Hue[0.0], PointSize[0.015]}, 
  BaseStyle -> {FontFamily -> "Helvetica", FontSize -> 18}]

Schluß

Hier ist Schluß.