BMonSphere.jpg (365 × 356 pixels, taille du fichier : 10 kio, type MIME : image/jpeg)

Ce fichier et sa description proviennent de Wikimedia Commons.

Description Brownian Motion on a Sphere. The generator of ths process is ½ times the Laplace-Beltrami-Operator
Date été 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.
Auteur Thomas Steiner
Autorisation
(Réutilisation de ce fichier)
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.
 
Cette infographie JPG a été créée avec Blender.
w:fr:Creative Commons
paternité partage à l’identique
Vous êtes libre :
  • de partager – de copier, distribuer et transmettre cette œuvre
  • d’adapter – de modifier cette œuvre
Sous les conditions suivantes :
  • paternité – Vous devez donner les informations appropriées concernant l'auteur, fournir un lien vers la licence et indiquer si des modifications ont été faites. Vous pouvez faire cela par tout moyen raisonnable, mais en aucune façon suggérant que l’auteur vous soutient ou approuve l’utilisation que vous en faites.
  • partage à l’identique – Si vous modifiez, transformez ou vous basez sur cet élément, vous devez distribuer votre contribution sous une license identique ou compatible à celle de l’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)

Légendes

Ajoutez en une ligne la description de ce que représente ce fichier
Brownian Motion on a Sphere, as a process generated by the Laplace-Beltrami-Operator

Éléments décrits dans ce fichier

dépeint

image/jpeg

f51c8d9194ca77a5c2a7d77d21c292e592d5c6c5

10 693 octet

356 pixel

365 pixel

Historique du fichier

Cliquer sur une date et heure pour voir le fichier tel qu'il était à ce moment-là.

Date et heureVignetteDimensionsUtilisateurCommentaire
actuel22 décembre 2013 à 20:53Vignette pour la version du 22 décembre 2013 à 20:53365 × 356 (10 kio)Olli NiemitaloCropped (in a JPEG-lossless way)
28 septembre 2007 à 23:53Vignette pour la version du 28 septembre 2007 à 23:53783 × 588 (14 kio)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

Usage global du fichier

Les autres wikis suivants utilisent ce fichier :