File:BMonSphere.jpg

From Wikipedia the free encyclopedia

BMonSphere.jpg(365 × 356 pixels, file size: 10 KB, MIME type: image/jpeg)

Description Brownian Motion on a Sphere. The generator of ths process is ½ times the Laplace-Beltrami-Operator
Date Summer 2007
date QS:P,+2007-00-00T00:00:00Z/9,P4241,Q40720564
(blender file as of 28.06.2007)
Source read some papers (eg Price, Gareth C.; Williams, David: "Rolling with “slipping”" : I. Séminaire de probabilités de Strasbourg, 17 (1983), p. 194-197 You can download it from http://www.numdam.org/item?id=SPS_1983__17__194_0) use the GNU R code and the python code (in blender3d) to create this image.
Author Thomas Steiner
Permission
(Reusing this file)
Thomas Steiner put it under the CC-by-SA 2.5. If you use the python code or the R code, please give a reference to Christian Bayer and Thomas Steiner.
 
This image was created with Blender.
w:en:Creative Commons
attribution share alike
This file is licensed under the Creative Commons Attribution-Share Alike 2.5 Generic license.
You are free:
  • to share – to copy, distribute and transmit the work
  • to remix – to adapt the work
Under the following conditions:
  • attribution – You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
  • share alike – If you remix, transform, or build upon the material, you must distribute your contributions under the same or compatible license as the original.

code

Perhaps you grab the source from the "edit" page without the wikiformating.

GNU R

This creates the paths and saves them into textfiles that can be read by blender. There are also paths for BMs on a torus.

# calculate a Brownian motion on the sphere; the output is a list # consisting of: # Z ... BM on the sphere # Y ... tangential BM, see Price&Williams # b ... independent 1D BM (see Price & Williams) # B ... generating 3D BM # n ... number of time-steps in the discretization # T ... the above processes are given on a uniform mesh of size #       n on [0,T] euler = function(x0, T, n) {   # initialize objects   dt = T/(n-1);   dB = matrix(rep(0,3*(n-1)),ncol=3, nrow=n-1);   dB[,1] = rnorm(n-1, 0, sqrt(dt));   dB[,2] = rnorm(n-1, 0, sqrt(dt));   dB[,3] = rnorm(n-1, 0, sqrt(dt));   Z = matrix(rep(0,3*n), ncol=3, nrow=n);   dZ = matrix(rep(0,3*(n-1)), ncol=3, nrow=n-1);   Y = matrix(rep(0,3*n), ncol=3, nrow=n);   B = matrix(rep(0,3*n), ncol=3, nrow=n);   b = rep(0, n);   Z[1,] = x0;    #do the computation   for(k in 2:n){     B[k,] = B[k-1,] + dB[k-1,];     dZ[k-1,] = cross(Z[k-1,],dB[k-1,]) - Z[k-1,]*dt;     Z[k,] = Z[k-1,] + dZ[k-1,];     Y[k,] = Y[k-1,] - cross(Z[k-1,],dZ[k-1,]);     b[k] = b[k-1] + dot(Z[k-1,],dB[k-1,]);   }   return(list(Z = Z, Y = Y, b = b, B = B, n = n, T = T)); }  # write the output from euler in csv-files euler.write = function(bms, files=c("Z.csv","Y.csv","b.csv","B.csv"),steps=bms$n){   bigsteps=round(seq(1,bms$n,length=steps))   write.table(bms$Z[bigsteps,],file=files[1],col.names=F,row.names=F,sep=",",dec=".");   write.table(bms$Y[bigsteps,],file=files[2],col.names=F,row.names=F,sep=",",dec=".");   write.table(bms$b[bigsteps],file=files[3],col.names=F,row.names=F,sep=",",dec=".");   write.table(bms$B[bigsteps,],file=files[4],col.names=F,row.names=F,sep=",",dec="."); }  # calculate a Brownian motion on a 3-d torus with outer # radius R and inner radius r eulerTorus = function(x0, r, R, t, n) {   # initialize objects   dt = t/(n-1);   dB = matrix(rep(0,3*(n-1)),ncol=3, nrow=n-1);   dB[,1] = rnorm(n-1, 0, sqrt(dt));   dB[,2] = rnorm(n-1, 0, sqrt(dt));   dB[,3] = rnorm(n-1, 0, sqrt(dt));   Z = matrix(rep(0,3*n), ncol=3, nrow=n);   B = matrix(rep(0,3*n), ncol=3, nrow=n);   dZ = matrix(rep(0,3*(n-1)), ncol=3, nrow=n-1);   Z[1,] = x0;   nT = rep(0,3);    #do the computation   for(k in 2:n){     B[k,] = B[k-1,] + dB[k-1,];     nT = nTorus(Z[k-1,],r,R);     dZ[k-1,] = cross(nT, dB[k-1,]) + HTorus(Z[k-1,],r,R)*nT*dt;     Z[k,] = Z[k-1,] + dZ[k-1,];   }   return(list(Z = Z, B = B, n = n, t = t)); }  # write the output from euler in csv-files torus.write = function(bmt, files=c("tZ.csv","tB.csv"),steps=bmt$n){   bigsteps=round(seq(1,bmt$n,length=steps))   write.table(bmt$Z[bigsteps,],file=files[1],col.names=F,row.names=F,sep=",",dec=".");   write.table(bmt$B[bigsteps,],file=files[2],col.names=F,row.names=F,sep=",",dec="."); }  # "defining" function of a torus fTorus = function(x,r,R){   return((x[1]^2+x[2]^2+x[3]^2+R^2-r^2)^2 - 4*R^2*(x[1]^2+x[2]^2)); }  # normal vector of a 3-d torus with outer radius R and inner radius r nTorus = function(x, r, R) {   c1 = x[1]*(x[1]^2+x[2]^2+x[3]^2-R^2-r^2)/(3*x[1]^4*x[2]^2+3*x[3]^4*x[2]^2     +3*x[3]^4*x[1]^2+6*x[3]^2*x[1]^2*x[2]^2+3*x[1]^2*x[2]^4+3*x[3]^2*x[2]^4     -2*x[3]^2*R^2*r^2-4*x[1]^2*x[2]^2*R^2+x[1]^6+x[2]^6+x[3]^6+3*x[3]^2*x[1]^4     -4*x[1]^2*x[2]^2*r^2-4*x[1]^2*x[3]^2*r^2+2*R^2*x[1]^2*r^2     -4*x[2]^2*x[3]^2*r^2+2*R^2*x[2]^2*r^2-2*x[1]^4*R^2-2*x[1]^4*r^2     +R^4*x[1]^2+x[1]^2*r^4-2*x[2]^4*R^2-2*x[2]^4*r^2+R^4*x[2]^2+x[2]^2*r^4     +x[3]^2*R^4+x[3]^2*r^4-2*x[3]^4*r^2+2*x[3]^4*R^2)^(1/2);   c2 = x[2]*(x[1]^2+x[2]^2+x[3]^2-R^2-r^2)/(3*x[1]^4*x[2]^2+3*x[3]^4*x[2]^2     +3*x[3]^4*x[1]^2+6*x[3]^2*x[1]^2*x[2]^2+3*x[1]^2*x[2]^4+3*x[3]^2*x[2]^4     -2*x[3]^2*R^2*r^2-4*x[1]^2*x[2]^2*R^2+x[1]^6+x[2]^6+x[3]^6     +3*x[3]^2*x[1]^4-4*x[1]^2*x[2]^2*r^2-4*x[1]^2*x[3]^2*r^2+2*R^2*x[1]^2*r^2     -4*x[2]^2*x[3]^2*r^2+2*R^2*x[2]^2*r^2-2*x[1]^4*R^2-2*x[1]^4*r^2+R^4*x[1]^2     +x[1]^2*r^4-2*x[2]^4*R^2-2*x[2]^4*r^2+R^4*x[2]^2+x[2]^2*r^4+x[3]^2*R^4     +x[3]^2*r^4-2*x[3]^4*r^2+2*x[3]^4*R^2)^(1/2);   c3 = (x[1]^2+x[2]^2+x[3]^2+R^2-r^2)*x[3]/(3*x[1]^4*x[2]^2+3*x[3]^4*x[2]^2                                             +3*x[3]^4*x[1]^2                                             +6*x[3]^2*x[1]^2*x[2]^2                                             +3*x[1]^2*x[2]^4+3*x[3]^2*x[2]^4                                             -2*x[3]^2*R^2*r^2                                             -4*x[1]^2*x[2]^2*R^2+x[1]^6                                             +x[2]^6+x[3]^6+3*x[3]^2*x[1]^4                                             -4*x[1]^2*x[2]^2*r^2                                             -4*x[1]^2*x[3]^2*r^2                                             +2*R^2*x[1]^2*r^2                                             -4*x[2]^2*x[3]^2*r^2                                             +2*R^2*x[2]^2*r^2-2*x[1]^4*R^2                                             -2*x[1]^4*r^2+R^4*x[1]^2                                             +x[1]^2*r^4-2*x[2]^4*R^2                                             -2*x[2]^4*r^2+R^4*x[2]^2                                             +x[2]^2*r^4+x[3]^2*R^4                                             +x[3]^2*r^4-2*x[3]^4*r^2                                             +2*x[3]^4*R^2)^(1/2);   return(c(c1,c2,c3)); }  # mean curvature of a 3-d torus with outer radius R and inner radius r HTorus = function(x, r, R){   return( -(3*x[1]^4*r^4+4*x[2]^6*x[3]^2+4*x[1]^6*x[2]^2-3*x[2]^4*x[3]^2*R^2             -2*x[1]^6*R^2+4*x[1]^2*x[3]^6+x[3]^6*R^2+4*x[2]^4*R^2*r^2-x[1]^2*r^6             -x[2]^2*r^6+x[2]^4*R^4+4*x[2]^2*x[3]^2*R^4+6*x[2]^2*x[3]^2*r^4             -2*x[1]^2*R^2*r^4-x[1]^2*R^4*r^2-9*x[1]^4*x[2]^2*r^2             -9*x[1]^4*x[3]^2*r^2+4*x[1]^4*R^2*r^2+12*x[1]^2*x[3]^4*x[2]^2             -3*x[2]^6*r^2+4*x[1]^6*x[3]^2+3*x[3]^4*r^4-x[3]^4*R^4             -9*x[2]^4*x[3]^2*r^2+2*x[2]^2*x[3]^2*R^2*r^2+4*x[1]^2*x[2]^6             -6*x[1]^2*x[3]^2*x[2]^2*R^2-x[3]^2*r^6+6*x[2]^4*x[3]^4+x[3]^8             +x[1]^8+x[2]^8-3*x[1]^6*r^2+6*x[1]^4*x[3]^4+12*x[1]^2*x[3]^2*x[2]^4             -6*x[1]^2*x[2]^4*R^2-2*x[3]^4*R^2*r^2-2*x[2]^2*R^2*r^4-x[2]^2*R^4*r^2             -9*x[2]^2*x[3]^4*r^2+x[3]^2*R^2*r^4+x[3]^2*R^4*r^2-9*x[1]^2*x[2]^4*r^2             +2*x[1]^2*R^4*x[2]^2+6*x[1]^2*x[2]^2*r^4-3*x[1]^4*x[3]^2*R^2             -6*x[1]^4*x[2]^2*R^2+4*x[1]^2*x[3]^2*R^4+6*x[1]^2*x[3]^2*r^4             -9*x[1]^2*x[3]^4*r^2+8*x[1]^2*R^2*x[2]^2*r^2+2*x[1]^2*x[3]^2*R^2*r^2             +x[1]^4*R^4-3*x[3]^6*r^2-2*x[2]^6*R^2+6*x[1]^4*x[2]^4-x[3]^2*R^6             -18*x[1]^2*x[2]^2*x[3]^2*r^2+4*x[2]^2*x[3]^6+12*x[1]^4*x[3]^2*x[2]^2             +3*x[2]^4*r^4)/(3*x[1]^4*x[2]^2+3*x[3]^4*x[2]^2+3*x[3]^4*x[1]^2                             +6*x[3]^2*x[1]^2*x[2]^2+3*x[1]^2*x[2]^4+3*x[3]^2*x[2]^4                             -2*x[3]^2*R^2*r^2-4*x[1]^2*x[2]^2*R^2+x[1]^6+x[2]^6                             +x[3]^6+3*x[3]^2*x[1]^4-4*x[1]^2*x[2]^2*r^2                             -4*x[1]^2*x[3]^2*r^2+2*R^2*x[1]^2*r^2                             -4*x[2]^2*x[3]^2*r^2+2*R^2*x[2]^2*r^2-2*x[1]^4*R^2                             -2*x[1]^4*r^2+R^4*x[1]^2+x[1]^2*r^4-2*x[2]^4*R^2                             -2*x[2]^4*r^2+R^4*x[2]^2+x[2]^2*r^4+x[3]^2*R^4                             +x[3]^2*r^4-2*x[3]^4*r^2+2*x[3]^4*R^2)^(3/2)); }  # calculate the cross product of the two 3-dim vectors # x and y. No argument-checking for performance reasons cross = function(x,y){   res = rep(0,3);   res[1] = x[2]*y[3] - x[3]*y[2];   res[2] = -x[1]*y[3] + x[3]*y[1];   res[3] = x[1]*y[2] - x[2]*y[1];   return(res); }  # calculate the inner product of two vectors of dim 3 # returns a number, not a 1x1-matrix! dot = function(x,y){   return(sum(x*y)); }  # calculate the cross product of the two 3-dim vectors # x and y. No argument-checking for performance reasons cross = function(x,y){   res = rep(0,3);   res[1] = x[2]*y[3] - x[3]*y[2];   res[2] = -x[1]*y[3] + x[3]*y[1];   res[3] = x[1]*y[2] - x[2]*y[1];   return(res); }  ############# ### main-teil set.seed(280180) et=eulerTorus(c(3,0,0),3,5,19,10000) torus.write(et,steps=9000) # #bms=euler(c(1,0,0),4,70000) #euler.write(bms,steps=10000)  

blender3d

The blender (python) code to create a image that looks almost like this one. Play around...

## import data from matlab-text-file and draw BM on the S^2  ## (c) 2007 by Christan Bayer and Thomas Steiner  from Blender import Curve, Object, Scene, Window, BezTriple, Mesh, Material, Camera, World from math import *  ##import der BM auf der Kugel aus einem csv-file def importcurve(inpath="Z.csv"):         infile = open(inpath,'r')         lines = infile.readlines()         vec=[]         for i in lines:                 li=i.split(',')                 vec.append([float(li[0]),float(li[1]),float(li[2].strip())])         infile.close()         return(vec)  ##function um aus einem vektor (mit den x,y,z Koordinaten) eine Kurve zu machen def vec2Cur(curPts,name="BMonSphere"):         bztr=[]         bztr.append(BezTriple.New(curPts[0]))         bztr[0].handleTypes=(BezTriple.HandleTypes.VECT,BezTriple.HandleTypes.VECT)         cur=Curve.New(name) ##TODO wenn es das Objekt schon gibt, dann nicht neu erzeugen         cur.appendNurb(bztr[0])         for i in range(1,len(curPts)):                 bztr.append(BezTriple.New(curPts[i]))                 bztr[i].handleTypes=(BezTriple.HandleTypes.VECT,BezTriple.HandleTypes.VECT)                 cur[0].append(bztr[i])         return( cur )  #erzeugt einen kreis, der später die BM umgibt (liegt in y-z-Ebene) def circle(r,name="tubus"):         bzcir=[]         bzcir.append(BezTriple.New(0.,-r,-4./3.*r, 0.,-r,0., 0.,-r,4./3.*r))         bzcir[0].handleTypes=(BezTriple.HandleTypes.FREE,BezTriple.HandleTypes.FREE)         cur=Curve.New(name) ##TODO wenn es das Objekt schon gibt, dann nicht neu erzeugen         cur.appendNurb(bzcir[0])         #jetzt alle weietren pkte         bzcir.append(BezTriple.New(0.,r,4./3.*r, 0.,r,0., 0.,r,-4./3.*r))         bzcir[1].handleTypes=(BezTriple.HandleTypes.FREE,BezTriple.HandleTypes.FREE)         cur[0].append(bzcir[1])         bzcir.append(BezTriple.New(0.,-r,-4./3.*r, 0.,-r,0., 0.,-r,4./3.*r))         bzcir[2].handleTypes=(BezTriple.HandleTypes.FREE,BezTriple.HandleTypes.FREE)         cur[0].append(bzcir[2])         return ( cur )  #erzeuge mit skript eine (glas)kugel (UVSphere) def sphGlass(r=1.0,name="Glaskugel",n=40,smooth=0):         glass=Mesh.New(name) ##TODO wenn es das Objekt schon gibt, dann nicht neu erzeugen         for i in range(0,n):                 for j in range(0,n):                         x=sin(j*pi*2.0/(n-1))*cos(-pi/2.0+i*pi/(n-1))*1.0*r                         y=cos(j*pi*2.0/(n-1))*(cos(-pi/2.0+i*pi/(n-1)))*1.0*r                         z=sin(-pi/2.0+i*pi/(n-1))*1.0*r                         glass.verts.extend(x,y,z)         for i in range(0,n-1):                  for j in range(0,n-1):                         glass.faces.extend([i*n+j,i*n+j+1,(i+1)*n+j+1,(i+1)*n+j])                         glass.faces[i*(n-1)+j].smooth=1         return( glass )  def torus(r=0.3,R=1.4):          krGro=circle(r=R,name="grTorusKreis")           #jetzt das material ändern def verglasen(mesh):         matGlass = Material.New("glas") ##TODO wenn es das Objekt schon gibt, dann nicht neu erzeugen         #matGlass.setSpecShader(0.6)         matGlass.setHardness(30) #für spec: 30         matGlass.setRayMirr(0.15)         matGlass.setFresnelMirr(4.9)         matGlass.setFresnelMirrFac(1.8)         matGlass.setIOR(1.52)         matGlass.setFresnelTrans(3.9)         matGlass.setSpecTransp(2.7)         #glass.materials.setSpecTransp(1.0)         matGlass.rgbCol = [0.66, 0.81, 0.85]         matGlass.mode |= Material.Modes.ZTRANSP         matGlass.mode |= Material.Modes.RAYTRANSP         #matGlass.mode |= Material.Modes.RAYMIRROR         mesh.materials=[matGlass]         return ( mesh )  def maleBM(mesh):         matDraht = Material.New("roterDraht") ##TODO wenn es das Objekt schon gibt, dann nicht neu erzeugen         matDraht.rgbCol = [1.0, 0.1, 0.1]         mesh.materials=[matDraht]         return( mesh )  #eine solide Mesh-Ebene (Quader) # auf der höhe ebh, dicke d, seitenlänge (quadratisch) 2*gr def ebene(ebh=-2.5,d=0.1,gr=6.0,name="Schattenebene"):         quader=Mesh.New(name) ##TODO wenn es das Objekt schon gibt, dann nicht neu erzeugen         #obere ebene         quader.verts.extend(gr,gr,ebh)         quader.verts.extend(-gr,gr,ebh)         quader.verts.extend(-gr,-gr,ebh)         quader.verts.extend(gr,-gr,ebh)         #untere ebene         quader.verts.extend(gr,gr,ebh-d)         quader.verts.extend(-gr,gr,ebh-d)         quader.verts.extend(-gr,-gr,ebh-d)         quader.verts.extend(gr,-gr,ebh-d)         quader.faces.extend([0,1,2,3])         quader.faces.extend([0,4,5,1])         quader.faces.extend([1,5,6,2])         quader.faces.extend([2,6,7,3])         quader.faces.extend([3,7,4,0])         quader.faces.extend([4,7,6,5])         #die ebene einfärben         matEb = Material.New("ebenen_material") ##TODO wenn es das Objekt schon gibt, dann nicht neu erzeugen         matEb.rgbCol = [0.53, 0.51, 0.31]         matEb.mode |= Material.Modes.TRANSPSHADOW         matEb.mode |= Material.Modes.ZTRANSP         quader.materials=[matEb]         return (quader)  ################### #### main-teil ####  # wechsel in den edit-mode editmode = Window.EditMode() if editmode: Window.EditMode(0)  dataBMS=importcurve("C:/Dokumente und Einstellungen/thire/Desktop/bmsphere/Z.csv") #dataBMS=importcurve("H:\MyDocs\sphere\Z.csv") BMScur=vec2Cur(dataBMS,"BMname") #dataStereo=importcurve("H:\MyDocs\sphere\stZ.csv") #stereoCur=vec2Cur(dataStereo,"SterName")  cir=circle(r=0.01)  glass=sphGlass() glass=verglasen(glass) ebe=ebene()  #jetzt alles hinzufügen scn=Scene.GetCurrent() obBMScur=scn.objects.new(BMScur,"BMonSphere") obcir=scn.objects.new(cir,"round") obgla=scn.objects.new(glass,"Glaskugel") obebe=scn.objects.new(ebe,"Ebene") #obStereo=scn.objects.new(stereoCur,"StereoCurObj")  BMScur.setBevOb(obcir) BMScur.update() BMScur=maleBM(BMScur)  #stereoCur.setBevOb(obcir) #stereoCur.update()  cam = Object.Get("Camera")  #cam.setLocation(-5., 5.5, 2.9)  #cam.setEuler(62.0,-1.,222.6) #alternativ, besser?? cam.setLocation(-3.3, 8.4, 1.7)  cam.setEuler(74,0,200)  world=World.GetCurrent() world.setZen([0.81,0.82,0.61]) world.setHor([0.77,0.85,0.66])  if editmode: Window.EditMode(1)  # optional, zurück n den letzten modus           #ergebnis von #set.seed(24112000) #sbm=euler(c(0,0,-1),T=1.5,n=5000) #euler.write(sbm)  

Captions

Brownian Motion on a Sphere, as a process generated by the Laplace-Beltrami-Operator

Items portrayed in this file

depicts

image/jpeg

f51c8d9194ca77a5c2a7d77d21c292e592d5c6c5

10,693 byte

356 pixel

365 pixel

File history

Click on a date/time to view the file as it appeared at that time.

Date/TimeThumbnailDimensionsUserComment
current19:53, 22 December 2013Thumbnail for version as of 19:53, 22 December 2013365 × 356 (10 KB)Olli NiemitaloCropped (in a JPEG-lossless way)
22:53, 28 September 2007Thumbnail for version as of 22:53, 28 September 2007783 × 588 (14 KB)Thire{{Information |Description = Brownian Motion on a Sphere |Source = read some papere ;) use the GNU R code and the python code (in blender3d) to create this image. |Date = summer 2007 (blender file as of ) |Author = Thomas Steiner |P
The following pages on the English Wikipedia use this file (pages on other projects are not listed):

Global file usage

The following other wikis use this file: