File:BMonSphere.jpg
 | This image was uploaded in the JPEG format even though it consists of non-photographic data. This information could be stored more efficiently or accurately in the PNG or SVG format. If possible, please upload a PNG or SVG version of this image without compression artifacts, derived from a non-JPEG source (or with existing artifacts removed). After doing so, please tag the JPEG version with {{Superseded|NewImage.ext}} and remove this tag. This tag should not be applied to photographs or scans. If this image is a diagram or other image suitable for vectorisation, please tag this image with {{Convert to SVG}} instead of {{BadJPEG}}. If not suitable for vectorisation, use {{Convert to PNG}}. For more information, see {{BadJPEG}}. |  |
| Description | Brownian Motion on a Sphere. The generator of ths process is ½ times the Laplace-Beltrami-Operator |
| Date | Summer 2007 (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 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)
Category:Images with R source code
Category:Brownian motion
Category:3D plots
Category:Single spheres (raster graphics)
Category:Images with Python source code
Category:Yellow-green background
Category:Wiener process
Category:3D plots
Category:Brownian motion
Category:CC-BY-SA-2.5
Category:Created with Blender
Category:Images with Python source code
Category:Images with R source code
Category:Images with inappropriate JPEG compression
Category:Single spheres (raster graphics)
Category:Wiener process
Category:Yellow-green background