Till Eulenspiegels lustige Serie/ Schwingende Platten/ Skripte
Kurzbeschreibung
BearbeitenDie Skripte sind auf Systemen wie dem Raspberry Pi lauffähig und sind autonom. Sie benutzen absolut keine speziellen Bibliotheken. Die Listings verzichten auf die farbenfrohe Auszeichnung von Syntax, damit man sie einfach aus der Webseite ausschneiden kann.
Das Skript ritz.py berechnet mit den Formeln von Ritz (sollten alle fehlerfrei abgetippt sein?) die Serien von Chladni-Klangfiguren als Svg-Dateien. Die lange Laufzeit von 1 Min wird in einem behinderten Algorithmus verbraucht, der Punktmengen zu Linien vereinen will und das nur schlecht hinkriegt. Lohnt sich kaum, den zu verbessern. Es entstehen 'single-xy.svg' für Einzelfiguren, sowie dekorativ bunte aber sonst nutzlose Ueberlagerungen 'multi-x.svg'. Massen von unverständlichen Debug-Zeilen werden aufs Terminal gespuckt. Alles ignorieren.
Angebaut wurde eine vollständige Neuberechnung der stehenden Wellen von Grund auf. Der Grad der Approximation wurde um Eins angehoben, denn wir haben Mahlwerke für Zahlen, von denen ein Autor vor 100 Jahren nur träumen konnte. Die Matrixelemente des Potenzial-Operators werden zwischen orthonormalen Sätzen von Basisfunktionen errechnet und die Eigenwerte/Eigenvektoren der symmetrischen Matrizen komplett rausgeholt. Der Algorithmus der Jacobi-Rotationen tut es blitzartig und genau; die Dimension der Matrizen bleibt unter 16. Es gibt mehr Graphen und deswegen 2 Minuten Laufzeit. Die Ergebnis-Datei "ritzmodes.txt" enthält alle Basisfamilien, Eigenwerte, Eigenvektoren. Nur wenige Graphen haben subtile Unterschiede zwischen beiden Versionen.
Das Skript chladni.py rechnet eine Approximation der Klangfiguren mit einem anderen Orthonormalsystem durch, nämlich mit den Polynomen an Stelle der physikalisch besser motivierten Balkenfunktionen. Aber mit einer etwas höheren Dimension der Matrizen kommt man auch hier zu ganz brauchbaren Ergebnissen. Die Integralrechnung für die Matrixelemente ist leichter. Chladni.py importiert Funktionen aus Ritz.py.
Un die Versionen im Vergleich anzuschauen, gibt es dieses Shellskript, das die Bilder paarweise verschmilzt und anzeigt. Man benutze Alt-F4, um weiter zu gehen.
###################
# watch.sh compare image pairs from 2 modes of ritz.py
# depends on: rsvg-convert pngtopnm pnmcat gview.
goon() {
local A F G L N
L=A; N=1
if [ -n "$1" ]; then L=$1; fi
while [ $N -lt 16 ]; do
A="$L$N"; F="single-1-$A.svg"; G="single-2-$A.svg"; N=$((N+1))
if [ -f $F ]; then
rsvg-convert $F -z 0.8 | pngtopnm >f.pnm
rsvg-convert $G -z 0.8 | pngtopnm >g.pnm
pnmcat -leftright f.pnm g.pnm >h.pnm
echo "Chladni $A"
gview h.pnm # ersetze durch xview,... irgendeinen Bildbetrachter
else N=20; fi
done
}
for X in A B C D E F; do goon $X ; done
##############
Listing ritz.py
Bearbeiten# -*- coding: utf-8 -*-
# W.Ritz 1909, Loesungen der Plattenschwingung. Stehende Wellen, Knotenlinien.
# quadratische Platte, Basis aus 7 mal 7 Wellenfunktionen
from math import sqrt, pi, sin, cos, sinh, cosh
def div(p,q) : return int(p/q)
def tanh(x): return sinh(x)/cosh(x)
def coth(x) : return cosh(x)/sinh(x)
###### primitive plot
class svgdump :
def plotbegin(sd, w,h,f) : # args: width height zoomfactor
W=str(w); H=str(h); u=''; fw=str(f*w); fh=str(f*h)
s=( '<?xml version="1.0" encoding="UTF-8" standalone="no"?>\n'+
'<!DOCTYPE svg>\n'+
'<svg xmlns="http://www.w3.org/2000/svg"\n'+
' xmlns:xlink="http://www.w3.org/1999/xlink"\n'+
' width="'+W+'" height="'+H+'" viewBox="0 0 '+fw+' '+fh+'">\n')
return s
def plotend(sd) : return '</svg>\n'
def move(sd,x,y) : return 'M'+str(x)+','+str(y)
def line(sd,x,y) : return 'L'+str(x)+','+str(y)
def ball(sd, x,y,r,col) :
return ('<circle fill="'+col+'" cx="'+str(x)+'" cy="'+str(y)+
'" r="'+str(r)+'" />\n')
def labl(sd,s,x,y, size=14,fill='#000') :
f=' font-family="Arial"'
t='<tspan x="'+str(x)+'" y="'+str(y)+'">'+s+'</tspan>'; fsz= str(size)
return ('<text fill="'+fill+'"'+f+' font-size="'+fsz+'px">'+t+'</text>\n')
def path(sd,w,s,lc='#000',fc='none') : # path w=with lc=linecolor fc=fillcol
t='<path fill="'+fc+'" stroke="'+lc+'" stroke-width="'+str(w)+'" '
t+= ' stroke-linecap="round" '
return t + 'd="\n'+s+'" />\n'
# end class svgdump
def aprx(x) : # 3-digit approximation
y= (-x) if(x<0) else x; z= int(1e3*y+1e-3); z= (-z) if (x<0) else z
return z/1000.0
def curve(sd,scale,xy, fill='#000', width=1, bkgr='none') :
(xs,fx,xmin, ys,fy,ymin) = scale; s=''; n= div(len(xy),2)
for i in range(n) :
x= aprx(xmin+fx*(xy[2*i]-xs)); y= aprx(ymin+fy*(xy[2*i+1]-ys))
s+= sd.move(x,y) if(i==0) else sd.line(x,y)
return sd.path(width,s,lc=fill,fc=bkgr)
def curveset(sd,scale,xyz, fill='#000', width=1, bkgr='none') :
(xs,fx,xmin, ys,fy,ymin) = scale; s=''; n= len(xyz)
print('curveset n='+str(n))
for i in range(n) :
px,py,pz= tuple(xyz[i])
x= aprx(xmin+fx*(px-xs)); y= aprx(ymin+fy*(py-ys))
s+= sd.move(x,y) if(pz==0) else sd.line(x,y)
return sd.path(width,s,lc=fill,fc=bkgr)
def ballset(sd,scale,plt, fill='#000', r=5) :
(xs,fx,xmin, ys,fy,ymin) = scale; s=''; n= len(plt)
for i in range(n) :
x= aprx(xmin+fx*(plt[i][0]-xs)); y= aprx(ymin+fy*(plt[i][1]-ys))
s+= sd.ball(x,y,r,fill)
return s
def label(sd,scale,txt,tx,ty) : # tx,ty plot coordinates not pixels
(xs,fx,xmin, ys,fy,ymin) = scale
x= aprx(xmin+fx*(tx-xs)); y= aprx(ymin+fy*(ty-ys))
return sd.labl(txt,x,y)
def bareplot(fn,x,y) : # common x list, numerous y lists. write file fn.svg
n=len(x); m=len(y); xys=[[]]*(m+1); winx=800; winy=600
xmin=min(x); xmax=max(x); ymin=y[0][0]; ymax=ymin
for k in range(m): z=min(y[k]); ymin=min(ymin,z)
for k in range(m): z=max(y[k]); ymax=max(ymax,z)
xyf=[ xmin,ymin, xmax,ymin, xmax,ymax, xmin,ymax, xmin,ymin]; xys[0]=xyf
for k in range(m) :
cv=[0]*(2*n); xys[k+1]= cv
for i in range(n) : cv[2*i]=x[i]; cv[2*i+1]=y[k][i]
xs,fx,xpix= xmin, (winx-50)/(xmax-xmin+0.0),25 # scale of user units
ys,fy,ypix= ymin,(-winy+50)/(ymax-ymin+0.0),winy-25
scale=(xs,fx,xpix, ys,fy,ypix)
colo=['#777','#000','#00a','#077','#0a0','#770','#a00']
sd=svgdump(); buf=sd.plotbegin(winx,winy,1); i=0; dy=(ymax-ymin)/20.0
for t in xys : buf += curve(sd,scale, t, fill=colo[i%7]); i +=1
g=open(fn+'.svg','w'); g.write(buf+sd.plotend()); g.close()
#### end plot
def basisfunc(n,x,par) : # Ritz orthonormalbasis
u=[0.0]*n; k,a,b,q= tuple(par)
u[0]=1.0/sqrt(2); u[1]=x*sqrt(3/2.0); even=True
for i in range(2,n) :
if even : z=k[i]*x; u[i]= (a[i]*cos(z)+b[i]*cosh(z))/q[i]
else : z=k[i]*x; u[i]= (a[i]*sin(z)+b[i]*sinh(z))/q[i]
even= not even
return u # Tupel Basisfunktionen(x)
def initbasis(n=11,dbg=1) : # par=initbasis() Parameter der Basisfunktionen
def fehler(z,sign) : return sin(z)*cosh(z)+sign*cos(z)*sinh(z)
def interpol(xa,xb, ya,yb) : # interpolate x at y=0, given ya*yb < 0
return xa + ya*(xa-xb)/(yb-ya+0.0)
def findzero(z,sign) :
epsi=1e-10; dz=0.1
za=z-3*dz; ya=fehler(za,sign); bracket=False; ok= False; count=0
while (not bracket)and(count<25) : # Klammer um Nullstelle
zb=za+dz; yb=fehler(zb,sign); bracket=(yb*ya)<0
if not bracket : za=zb;ya=yb; count+=1
if bracket :
emax=epsi*(abs(ya)+abs(yb)); # print('bracket '+str([count,ya,yb]))
while (not ok) and (count<50) :
z=interpol(za,zb,ya,yb); y=fehler(z,sign)
if (ya*y)>0 : ya=y;za=z
else : yb=y;zb=z
ok= (abs(y)<emax); count +=1
if ok and (dbg==1): print('n,c,z='+str([n,count,z]))
if not ok: print('exit not ok n='+str(n)+' count='+str(count)); exit()
return z
###
k=[0]*n; a=[0]*n; b=[0]*n; q=[0]*n; even=True
for m in range(2,n) :
sign= 1 if even else -1; z= (m-1/2.0)*pi/2.0
z= findzero(z,sign) # solve tan(z)+ sign*tanh(z)=0
if even : a[m]=cosh(z); b[m]=cos(z)
else : a[m]=sinh(z); b[m]=sin(z)
k[m]=z; q[m]=sqrt(a[m]*a[m]+b[m]*b[m])
even= not even
return [k,a,b,q]
def platewave(tpe,x,y,pars) : # alle Ritz-Approximationen der stehenden Wellen
u=basisfunc(7,x,pars); u0,u1,u2,u3,u4,u5,u6= tuple(u)
v=basisfunc(7,y,pars); v0,v1,v2,v3,v4,v5,v6= tuple(v)
# print('u='+str(u)); print('v='+str(v))
deltamu=0.0; lambd=[0]*15; w=[0]*15; wb=[0]*15
# A. in x und y ungerade, symmetrisch: w(x,y)=w(y,x)=-w(-x,y)=-w(x,-y)
if tpe=='A' :
lambd[0]= 12.43 - 18.0*deltamu
w[0]= (u1*v1 + 0.0394*(u1*v3+v1*u3)-0.0040*u3*v3-0.0034*(u1*v5+u5*v1)
+0.0011*(u3*v5+u5*v3)-0.0019*u5*v5)
lambd[1]=378-57*deltamu
w[1]= (-0.075*u1*v1+(u1*v3+u3*v1)+0.173*u3*v3+0.045*(u1*v5+u5*v1)
-0.015*(u3*v5+u5*v3)-0.029*u5*v5)
lambd[2]=1554
w[2]= (0.009*u1*v1-0.075*(u1*v3+v1*u3)+u3*v3-0.057*(u1*v5+u5*v1)
+0.121*(u3*v5+u5*v3)-0.007*u5*v5)
lambd[3]=2945
w[3]= u1*v5+u5*v1
lambd[4]=6303
w[4]= u3*v5+u5*v3
lambd[5]=13674
w[5]= u5*v5; last=5
# B. in x und y ungerade, antisymmetrisch: w(x,y)=-w(y,x)=-w(-x,y)=-w(x,-y)
elif tpe=='B' :
lambd[0]=316.1-270*deltamu
w[0]= u1*v3-v1*u3+ 0.0002*(u1*v5-v1*u5)+0.0033*(u3*v5-v3*u5)
lambd[1]=2713
w[1]=u1*v5-v1*u5
lambd[2]=5570
w[2]=u3*v5-v3*u5; last=2
# C. in x und y gerade, symmetrisch: w(x,y)=w(y,x)=w(-x,y)=w(x,-y)
elif tpe=='C' :
lambd[0]= 35.73+20.8*deltamu
w[0]=(u0*v2+u2*v0-0.0238*u2*v2+0.0130*(u0*v4+v0*u4)
+0.0026*(u2*v4+v2*u4)+0.0016*u4*v4)
lambd[1]=266.0-274*deltamu
w[1]= (0.0122*(u0*v2+v0*u2)+u2*v2-0.0188*(u0*v4+v0*u4)
+0.0880*(u4*v2+v4*u2)-0.0044*u4*v4)
lambd[2]=941
w[2]= u0*v4+v0*u4
lambd[3]=2020
w[3]= u2*v4+v2*u4
lambd[4]=5480
w[4]= u4*v4
lambd[5]=5640
w[5]= u0*v6+v0*u6
lambd[6]=7840
w[6]= u2*v6+v2*u6
lambd[7]=15120
w[7]= u4*v6+v4*u6
lambd[8]=28740
w[8]= u6*v6; last=8
# D. in x und y gerade, antisymm: w(x,y)=-w(y,x)=w(-x,y)=w(x,-y)
elif tpe=='D' :
lambd[0]=26.40
w[0]=u0*v2-v0*u2-0.0129*(u0*v4-v0*u4)-0.0045*(u2*v4-v2*u4)
lambd[1]= 886
w[1]=u0*v4-v0*u4
lambd[2]=1702
w[2]=u2*v4-v2*u4
lambd[3]=5500
w[3]= u0*v6-v0*u6
lambd[4]=7310
w[4]= u2*v6-v2*u6
lambd[5]=13840
w[5]= u4*v6-v4*u6; last=5
# E. Doppeltoene (Chladni-Figuren nicht 90Grd-symm,)
elif (tpe=='E')or(tpe=='F')or(tpe=='G') :
# Entartung a*w(x,y)+b*w(y,x) (a,b)=(1,0),(1,-1)
lambd[0]=80.8-73*deltamu
w[0]=(u1*v2-0.0682*u3*v0 +0.0760*u3*v2+0.0260*u1*v4
+0.0073*u5*v0-0.0027*u3*v4-0.0112*u5*v2+0.0030*u5*v4)
z= (v1*u2-0.0682*v3*u0 +0.0760*v3*u2+0.0260*v1*u4
+0.0073*v5*u0-0.0027*v3*u4-0.0112*v5*u2+0.0030*v5*u4)
wb[0]=w[0]-z
lambd[1]=237.1
w[1]= (+0.0678*u1*v2+u3*v0-0.0150*u3*v2+0.0355*u1*v4
+0.0000*u5*v0+0.0100*u3*v4-0.0007*u5*v2+0.0016*u5*v4)
z= (+0.0678*v1*u2+v3*u0-0.0150*v3*u2+0.0355*v1*u4
+0.0000*v5*u0+0.0100*v3*u4-0.0007*v5*u2+0.0016*v5*u4)
wb[1]=w[1]-z
lambd[2]=746
w[2]=(-0.0709*u1*v2+0.0214*u3*v0+u3*v2-0.1260*u1*v4
-0.0038*u5*v0+0.1234*u3*v4-0.0095*u5*v2-0.0100*u5*v4)
z= (-0.0709*v1*u2+0.0214*v3*u0+v3*u2-0.1260*v1*u4
-0.0038*v5*u0+0.1234*v3*u4-0.0095*v5*u2-0.0100*v5*u4)
wb[2]=w[2]-z
lambd[3]=1131
w[3]= u1*v4; wb[3]=u1*v4-u4*v1
lambd[4]=2497
w[4]=u5*v0; wb[4]=u5*v0-u0*v5
lambd[5]=3240
w[5]=u3*v4; wb[5]=u3*v4-u4*v3
lambd[6]=3927
w[6]=u5*v2; wb[6]=u5*v2-u2*v5
lambd[7]=9030
w[7]=u5*v4; wb[7]=u5*v4-v5*u4
lambd[8]=6036
w[8]=u1*v6; wb[8]=u1*v6-v1*u6
lambd[9]=10380
w[9]=u3*v6; wb[9]=u3*v6-v3*u6
lambd[10]=20400
w[10]=u5*v6; wb[10]=u5*v6-v5*u6; last=10
if (tpe=='F')or(tpe=='G') : w=wb
else : print('Invalid type tpe='+tpe); exit()
return lambd[:last+1], w[:last+1]
def ballplot(fn,xa,ya,xb,yb,pball,pcurve) : # dots in pball, lines in pcurve
winx=1000; winy=1000; xmin=xa;ymin=ya;xmax=xb;ymax=yb; rad=5; wid=3
xyf=[xmin,ymin, xmax,ymin, xmax,ymax, xmin,ymax, xmin,ymin]
xs,fx,xpix= xmin, (winx-50)/(xmax-xmin+0.0),25 # scale user units, pixel@min
ys,fy,ypix= ymin,(-winy+50)/(ymax-ymin+0.0),winy-25
scale=(xs,fx,xpix, ys,fy,ypix)
colo=['#777','#000','#00a','#077','#0a0','#770','#a00']
sd=svgdump(); buf=sd.plotbegin(winx,winy,1); i=0; dy=(ymax-ymin)/20.0
# for t in xys : buf += curve(sd,scale, t, fill=colo[i%7]); i +=1
buf += curve(sd,scale, xyf, fill=colo[i%7],bkgr='#ff9') # frame
for k in range(len(pball)) :
buf+= ballset(sd,scale,pball[k],fill=colo[(k+1)%7], r=rad)
for k in range(len(pcurve)) :
buf+= curveset(sd,scale,pcurve[k],fill=colo[(k+1)%7], width=wid)
g=open(fn+'.svg','w'); g.write(buf+sd.plotend()); g.close()
def getlines(plt,d,ang) : # transform pointset to lines, max neighb distance d
def dist(a,b,c,d) : u=c-a; v=d-b; return sqrt(u*u+v*v)
def inline(plt,i,k) : # forbid bend angles with cosine < ang.
j=prev[i]; ok= (j<0)
if not ok :
xa,ya=tuple(plt[j]); xb,yb=tuple(plt[i]); xc,yc=tuple(plt[k])
sp=(xc-xb)*(xb-xa)+(yc-yb)*(yb-ya) # scalar product
nab=dist(xa,ya, xb,yb); nbc=dist(xb,yb, xc,yc)
ok=(sp>(ang*nab*nbc))
return ok
n=len(plt); next=[-1]*n; prev=[-1]*n; seen=[0]*n; done=False; i=-1
print('Point set size '+str(n))
while not done :
if i<0 :
i=0 # find unseen isolated point
# while (i<n) and (seen[i]>0) and ((next[i]>=0)or(prev[i]>=0)) : i+=1
while (i<n) and ((seen[i]>0) or (next[i]>=0) or (prev[i]>=0)) : i+=1
# if i<n : seen[i]=1; # print('new start '+str(i)+'/'+str(n))
if i<n : # try to continue point i
xa,ya= tuple(plt[i]); k=0; dmin=1000; kmin=-1; seen[i]=1
while (k<n) and (next[i]<0) :
if (k!= i) and (prev[k]==-1) :
xb,yb= tuple(plt[k]); dx=dist(xa,ya,xb,yb)
if dx<dmin : dmin=dx; kmin=k
k+=1
if (dmin<d) and inline(plt,i,kmin) :
k=kmin; next[i]=k;prev[k]=i; ok = (xb>xa) or (yb>ya)
if (not ok) and (prev[i]<0) and(next[k]<0) : # flip orientation
next[k]=i;prev[i]=k
else : i=k
else : i=-1
else : done=True
single=0; loops=0; paths=[]; seen=[0]*n
for i in range(n) :
if (next[i]<0)and(prev[i]<0) : single+=1
elif seen[i]==0 :
seen[i]=1; links=0; j=i; go=True
while (next[j]>=0)and go : j=next[j]; go=(seen[j]==0); seen[j]=1;links+=1
k=i; go=True
while (prev[k]>=0)and go : k=prev[k]; go=(seen[k]==0); seen[k]=1;links+=1
print('Path '+str([k,j,links])); paths+=[[k,j]]
print('Singles: '+str(single)+' Paths: '+str(len(paths)))
return paths, next
def mirrors(xyz,subset) : # symmetrie +-x, +-y. fuer subsets ABCDE
# subset F,G: xyz doppelt so gross
n=len(xyz); f= 2 if(subset=='F') else 4; xyzm=[[]]*(f*n)
for i in range(n) :
x,y,z= tuple(xyz[i]); xyzm[i]=[x,y,z]; xyzm[n+i]=[-x,-y,z]
if subset != 'F' : xyzm[2*n+i]=[-x,y,z]; xyzm[3*n+i]=[x,-y,z]
return xyzm
def plotstuff(subset,curves,debug,version='1') :
if debug<=2 :
ballplot('multi-'+version+'-'+subset, -100,-100, 100,100, [], curves)
if debug==2 : # Einzelkurven
for i in range(len(curves)) :
ballplot('single-'+version+'-'+subset+str(i+1),
-100,-100, 100,100, [], [curves[i]])
def knotenbild(subset, w, debug=0, version='1') : # Chladni-Figur Knotenlinien
# debug=2 macht Kurven mit Namen 'single'+subset+version
dmin= 2.1; angle=0.6
n=len(w); m=len(w[0][0]); plot=[[]]*m; xy=[[]]*m
print('Knotenbild n='+str(n)+' m='+str(m)+' n2='+str(len(w[0])))
for k in range(m) : # Nullstellen in Teilmenge k
plt=[]
for i in range(n) :
z=w[0][i][k]
if z==0 : plt+=[[0,i]] # abs(z)<epsilon?
for i in range(1,n) :
z=w[i][0][k]
if z==0 : plt+= [[i,0]]
for j in range(1,n) : # lineare Interpolationen der Nulldurchgaenge
z=w[i][j][k]; x=w[i-1][j][k]; y=w[i][j-1][k]; u=w[i-1][j-1][k]
if z==0 : plt+= [[i,j]] # max 1 Punkt per Paar (i,j)
elif (z*x)<0 : plt += [[i-1 -x/(z-x+0.0), j]]
elif (z*y)<0 : plt += [[i,j-1 -y/(z-y+0.0)]]
# diagonal, z*u und x*y ? mehr Punkte
elif (u*z)<0 : plt += [[i-1 -u/(z-u+0.0), j-1 -u/(z-u+0.0) ]]
elif (x*y)<0 : plt += [[i-1 -x/(y-x+0.0), j-1 -y/(x-y+0.0) ]]
plot[k]=plt; print('k='+str(k)+' len(plt)='+str(len(plt)))
if debug==1 : ballplot('ball-'+subset,0,0,100,100,plot,[]) # 500K svg
curves=[]; # m= min(3,m) # debug
for i in range(m) :
paths,next = getlines(plot[i], dmin,angle) # schlecht und langsam
xyz=[]; plt=plot[i]
print('Plotting paths '+str(len(paths)))
for k in range(len(paths)) :
h,j=paths[k]; z=0
while h != j : xyz+=[[plt[h][0],plt[h][1],z]]; z=1; h=next[h]
xyz+=[[plt[j][0],plt[j][1],1]]
curves+= [xyz] if (debug==3) else [mirrors(xyz,subset)]
plotstuff(subset,curves,debug,version)
return curves
def simplecase(subset, dbg=0) : # subsets: A B C D E F G
# case F,G needs w on x=-100..100,y=0..100, later node[-x,-y]=node[x,y]
def vscale(f,v) : # scale vector v
n=len(v); u=[0]*n
for i in range(n) : u[i]=f*v[i]
return u
pars=initbasis(); nx=100; ny=100; w=[[]]*(nx+1)
nodiag=('EFG'.find(subset)>=0); xs= (-1) if(subset=='G') else 1
sign=(-1) if('BD'.find(subset)>=0) else 1 # EFG no x=y mirror!
for ix in range(nx+1) :
x= ix/(nx+0.0); w[ix]=[[]]*(ny+1); ymax= (nx+1) if nodiag else (ix+1)
for iy in range(ymax) :
y= iy/(ny+0.0)
la, w[ix][iy] = platewave(subset, xs*x,y,pars)
if not nodiag : # x=y mirror
for ix in range(nx+1) :
for iy in range(ix) : w[iy][ix]= vscale(sign,w[ix][iy]) # Symm-Achse x=y
print('Wellenwerte: '+str([len(w),len(w[0]),len(w[0][0])]))
if ('ABCDE'.find(subset)>=0) : curves=knotenbild(subset, w, debug=dbg)
if (subset=='F')or(subset=='G') : curves=knotenbild(subset, w, debug=3)
return curves
def specialcase(compute=0,version='1') :
def merger(c1,c2) :
n1=len(c1); n2=len(c2); d1=[[]]*(2*n1); d2=[[]]*(2*n2)
for i in range(n1) : x,y,z= tuple(c1[i]); d1[i]=[x,y,z]; d1[n1+i]=[-x,-y,z]
for i in range(n2) : x,y,z= tuple(c2[i]); d2[i]=[-x,y,z]; d2[n2+i]=[x,-y,z]
return d1+d2
if compute==0 :
curv1= simplecase('F') # quadrant x>0 y>0
curv2= simplecase('G') # quadrant x<0 y>0
else :
curv1= computedcase('F', big=(compute==2))
curv2= computedcase('G', big=(compute==2))
m1=len(curv1); m2=len(curv2); curves=[[]]*m1
if m1 != m2 : print('specialcase m1 != m2, exit'); exit()
for i in range(m1) : curves[i]= merger(curv1[i],curv2[i])
plotstuff('F', curves, debug=2, version=version)
### Eigenwerte nach wikipedias Algorithmus
def rotjacobi(a,x,count,dbg=0) : # after en.wikipedia, Jacobi_rotation
n=len(a); gain=0; tiny=1e-33; giant=1e33
for k in range(n): # upper triangle k<l, one similarity transf per pair
for l in range(k+1,n) : # todo: keep list of pivots= largest z per row?
z=a[k][l]; ok=(abs(z)<tiny)
if ok : a[k][l]=0 # experiment
else : # later transform with (c-s|sc) on indexpairs (kl)
gain+=abs(z) # sign t=-abs(b)+... !!
v=(a[l][l]-a[k][k]); b=v/(2.0*z); t= -abs(b)+sqrt(1+b*b)
t=(-t) if(b<0) else t # if v=0, t=1, angle pi/4.
if t>giant : s=0; c=1; tz=v
elif t<-giant : s=0; c=-1; tz=-v
else : u=sqrt(1.0+t*t); c=1.0/u; s=c*t; tz=t*z
a[k][k] -= tz; a[l][l] += tz; a[k][l]=0
for h in range(n) :
if (h!=k) and (h!=l) : # rotate pair {a_hk,a_hl}, use index symmetry.
(hk,kh)= (h,k) if (h<k) else (k,h); ak=a[hk][kh]
(hl,lh)= (h,l) if (h<l) else (l,h); al=a[hl][lh]
a[hk][kh]= c*ak-s*al; a[hl][lh]= s*ak+c*al
for h in range(n) : # rotate columns, no transpose later
xk=x[k][h];xl=x[l][h]; x[k][h]= c*xk-s*xl; x[l][h]= s*xk+c*xl
return gain
def eigenloop(ma) : # symmetric matrix ma
n=len(ma); x=[[]]*n; m=[[]]*n; epsi=1e-33; j=0; gain=1; eval=[0]*n
for i in range(n) :
x[i]=[0]*n; x[i][i]=1; m[i]=[0]*n
for k in range(n) : m[i][k]=ma[i][k]
while (j<(2*n))and(gain>epsi) : # matrix m is transformed
gain=rotjacobi(m,x,j,dbg=1); j+=1; # print('Gain='+str(gain))
for i in range(n) : eval[i]=m[i][i]
return eval,x # eval=eigenval list, rows of x are the eigenvectors
def eigentest(m,ev,x,dbg=0) : # Gleichungen m x_i = ev_i x_i
n=len(m); y=[0]*n; z=[0]*n; norms=0
for j in range(n) : # test eigenvector j, eigenval j
norm=0
for i in range(n) :
z[i]= ev[j]*x[j][i]; y[i]=0
for k in range(n) : y[i]+= m[i][k]*x[j][k]
d=y[i]-z[i]; norm += d*d
norms+= norm
if dbg==1 : print('Eigenval '+str(j)+' Norm='+str(norm))
return norms
#####
def fe(x) : return ('%11d' % x) if(type(x)==type(0)) else ('%+6.4e' % x)
def fdump(s,m) : # vector or matrix debug dump
n=len(m); ismat= (type(m[0])==type([])); t=''
if not ismat :
t='Vector '+s+' Dim='+str(n)+'\n'
for k in range(n) : t+=' '+fe(m[k])
else :
j=len(m[0]); t='Matrix '+s+' Dim='+str(n)+'*'+str(j)+'\n'
for i in range(n) :
for k in range(j) : t+=' '+fe(m[i][k])
if i<(n-1) : t+='\n'
return t
def norm(a) :
n=len(a); z=0.0
for i in range(n) : z += a[i]*a[i]
return sqrt(z)
def wsum(a,b,f,g) : # weighted sum of 2 vectors
n=len(a); c=[0]*n
for i in range(n) : c[i]= f*a[i]+g*b[i]
return c
def dbasisfunc(n,x,par,dg=1) : # Ritz orthonormalbasis, derivatives
du=[0.0]*n; k,a,b,q= tuple(par)
du[0]=0; du[1]=sqrt(3/2.0); even=True; y=1
if dg==1 : # first deriv
for i in range(2,n) :
if even : z=k[i]*x; du[i]= (-a[i]*sin(z)+b[i]*sinh(z))*k[i]/q[i]
else : z=k[i]*x; du[i]= ( a[i]*cos(z)+b[i]*cosh(z))*k[i]/q[i]
even= not even
else : # other degree up to 4
du[1]=0; feven=[cos,sin,cos,sin,cos]; fodd=[sin,cos,sin,cos,sin]
se=[1,-1,-1,1,1]; so=[1,1,-1,-1,1]
heven=[cosh,sinh,cosh,sinh,cosh]; hodd=[sinh,cosh,sinh,cosh,sinh]
for i in range(2,n) : # dg==0 function itself
y=1; z=k[i]*x
for j in range(dg) : y *= k[i]
if even : du[i]= (a[i]*se[dg]*feven[dg](z)+b[i]*heven[dg](z))*y/q[i]
else : du[i]= (a[i]*so[dg]*fodd[dg](z)+b[i]*hodd[dg](z))*y/q[i]
even= not even
return du # Tupel
def alphaomega(m,n,par) : # alpha= um'*un', omega=um"*un integrals.
if (m==0) or (n==0) : alpha=0; omega=0
elif (m==1)and(n==1) : alpha=3; omega=0
else :
k, u1, du1= tuple(par) # u1, du1; Liste von Basisfunc und Ableitg bei x=1
km=k[m]; kn=k[n]
if (m==n)and(m%2==0) :
ch=cosh(km); co= cos(km); ch2=ch*ch; co2=co*co; q=ch2+co2
alpha= (km*km*(ch2-co2)+6*km*(co2*ch2*tanh(km)))/q
omega= (-km*km*(ch2-co2)+2*km*(co2*ch2*tanh(km)))/q
elif (m==n) and (m%2==1) :
sh=sinh(km); si= sin(km); sh2=sh*sh; si2=si*si; q=sh2-si2
alpha= (km*km*(sh2+si2)+6*km*(si2*sh2*coth(km)))/q
omega= (-km*km*(sh2+si2)+2*km*(si2*sh2*coth(km)))/q
elif (m+n)%2 == 1 : alpha=0; omega=0 # equal parity --> vanishes
else : # vanish if m+n is odd
km4=km*km*km*km; kn4=kn*kn*kn*kn; q=km4-kn4
alpha= 2*(km4*u1[m]*du1[n]-kn4*u1[n]*du1[m])/q
omega= 2*km4*(du1[m]*u1[n]-u1[m]*du1[n])/q
return alpha,omega
def vmatrix(m,n,p,q, mu, par): # Element zwischen um*vn, up*vq.
if (m==p)and(n==q) :
k,u1,du1= tuple(par)
z=k[m];km4=z*z*z*z; z=k[n]; kn4=z*z*z*z
amm,omm= alphaomega(m,m,par); ann,onn= alphaomega(n,n,par)
v= km4+kn4+2*mu*omm*onn + 2*(1-mu)*amm*ann
else :
amp,omp= alphaomega(m,p,par); anq,onq= alphaomega(n,q,par)
aqn,oqn= alphaomega(q,n,par); apm,opm= alphaomega(p,m,par)
v= mu*(omp*oqn+opm*onq)+2*(1-mu)*anq*amp
return v
def ritzmatrix(serie='A', mode=0, normalize=False, big=False) :
# Auswahl einer symmetrischen Orthonormalbasis aus dem VONS der Stabwellen
# return: Matrix des V-operators in dieser Basis, Index-triplets der Basis
# dbg mode=1 : zeilen = Faktoren u1v1 u1v3 u1v5 u3v3 u3v5 u5v5 wie im Artikel
# big=False # big Approximation einen Grad hoeher
n=11; mu=0.225 # Querkontraktions-Parameter, Glas
pars=initbasis(n=n,dbg=0); k,a,b,q= tuple(pars)
u1= basisfunc(n,1,pars); du1=dbasisfunc(n,1,pars)
vpar=[k, u1, du1] # ixtrip= pairs with target index and sign contrib
if serie=='A' : # index triplets ui*vj itarget(base 1!). odd odd sym
ixtrip= [[1,1,1],[1,3,2],[3,1,2],[3,3,3],[1,5,4],[5,1,4],
[3,5,5],[5,3,5],[5,5,6]]
if big : ixtrip+=[[1,7,7],[7,1,7],[3,7,8],[7,3,8],[5,7,9],[7,5,9],[7,7,10]]
if serie=='B' : # odd odd antisym
ixtrip=[[1,3,1],[3,1,-1],[1,5,2],[5,1,-2],[3,5,3],[5,3,-3]]
if big: ixtrip+= [[1,7,4],[7,1,-4], [3,7,5],[7,3,-5], [5,7,6],[7,5,-6]]
if serie=='C' : # even even sym
ixtrip=[[0,2,1],[2,0,1],[2,2,2],[0,4,3],[4,0,3],[2,4,4],[4,2,4],[4,4,5],
[0,6,6],[6,0,6],[2,6,7],[6,2,7],[4,6,8],[6,4,8],[6,6,9]]
if big : ixtrip+=[[0,8,10],[8,0,10],[2,8,11],[8,2,11],[4,8,12],[8,4,12],
[6,8,13],[8,6,13],[8,8,14]]
if serie=='D' : # even even antisym
ixtrip=[[0,2,1],[2,0,-1],[0,4,2],[4,0,-2],[2,4,3],[4,2,-3],
[0,6,4],[6,0,-4],[2,6,5],[6,2,-5],[4,6,6],[6,4,-6]]
if big : ixtrip+= [[0,8,7],[8,0,-7],[2,8,8],[8,2,-8],[4,8,9],[8,4,-9],
[6,8,10],[8,6,-10]]
if serie=='E' : # odd even
ixtrip=[[1,2,1],[3,0,2],[3,2,3],[1,4,4],[5,0,5],[3,4,6],
[5,2,7],[5,4,8]]
if big : ixtrip+=[[1,6,9],[3,6,10],[5,6,11],[7,6,12],
[7,0,13],[7,2,14],[7,4,15]]
if (serie=='F')or(serie=='G') : # (odd even - even odd)
ixtrip=[[1,2,1],[2,1,-1], [3,0,2],[0,3,-2], [3,2,3],[2,3,-3],
[1,4,4],[4,1,-4], [5,0,5],[0,5,-5], [3,4,6],[4,3,-6],
[5,2,7],[2,5,-7], [5,4,8],[4,5,-8]]
if big : ixtrip+=[[1,6,9],[6,1,-9], [3,6,10],[6,3,-10], [5,6,11],[6,5,-11],
[7,6,12],[6,7,-12], [7,0,13],[0,7,-13],
[7,2,14],[2,7,-14], [7,4,15],[4,7,-15]]
rowtrip= ixtrip; coltrip=ixtrip # row and column triplets ui,vj,itarget
if (serie=='A')and(mode==1) :
rowtrip=[[1,1,1],[1,3,2],[3,3,3],[1,5,4],[3,5,5],[5,5,6]]
mcol=len(coltrip); mrow=len(rowtrip)
v=[[]]*mrow; mmat=0 # mmat size of final matrix
# ixpairs is finite-dim subspace in one of the symmetry classes
for i in range(mrow) :
v[i]=[0]*mcol; mi,ni,xi= tuple(rowtrip[i]); mmat=max(mmat,abs(xi))
for j in range(mcol) :
mj,nj,xj= tuple(coltrip[j]); v[i][j]= vmatrix(mi,ni,mj,nj, mu,vpar)
rm=[[]]*mmat; nm=[[]]*mmat # Ritz matrix and normalize count
for h in range(mmat) : rm[h]=[0]*mmat; nm[h]=[0]*mmat
for i in range(mrow) : # add to row i,xi
mi,ni,xi= tuple(rowtrip[i]); (si,xi) =(1,xi-1) if (xi>0) else (-1,-xi-1)
for j in range(mcol) : # addto column j,xj
mj,nj,xj= tuple(coltrip[j]); (sj,xj) =(1,xj-1) if (xj>0) else (-1,-xj-1)
rm[xi][xj]+= si*sj*v[i][j]; nm[xi][xj] +=1
if normalize : # divide rm by square-root of counts nm, for orthonormal
for i in range(mmat) :
for j in range(mmat) :
if nm[i][j]>1 : rm[i][j]= rm[i][j]/sqrt(nm[i][j])
return rm, ixtrip, pars # pars= Parameter des VONS
if __name__ == '__main__' : # do not export what follows
def eigenwave(ev,ixt,x,y,pars) : # stehende Wellen, vektoren ev in Basis ixt
u=basisfunc(9,x,pars); v=basisfunc(9,y,pars) # max Ordnung 8
bw=[0.0]*16; kmax=0; n=len(ixt)
for i in range(n) : # Basis-Wellenwerte
j,k,nd= tuple(ixt[i]); (sgn,ix)= (1,nd-1) if (nd>0) else (-1,-nd-1)
bw[ix]+= sgn*u[j]*v[k]; kmax=max(kmax,ix)
m=kmax+1; w=[0.0]*m # linearkombiniert zu Eigenwellen
for i in range(m) :
for j in range(m) : w[i]+= ev[i][j]*bw[j]
return w # m Ergebnisse am Punkt (x,y)
def dumpmodes(subset, ixt, ewert,evec,ref) : # Doku text output
n=len(ixt); m=len(evec); bw=['']*m; count=[0]*m
d='Serie '+subset+'. Basisfunktionen:'+str(m)+'. Faktor q= 1/sqrt(2)\n'
for i in range(n) : # Basis-Wellenformeln
j,k,nd= tuple(ixt[i]); (sgn,ix)= ('+',nd-1) if (nd>0) else ('-',-nd-1)
bw[ix]+= sgn+'u'+str(j)+'*v'+str(k); count[ix]+=1
for j in range(m) :
s=('('+bw[j]+')') if(count[j]==1) else ('('+bw[j]+')*q'); d+=(s+'\n')
d+='Eigenwerte und Eigenvektor-Koeffizienten in der Basis\n'
for j in range(m) :
s=fe(ewert[j])+' : ['
for k in range(m): s+=' '+ fe(evec[j][k])
d+=s+' ]\n'
d+= '**** Ritz Computer (Eigenwert-Vergleich)\n'; nr=len(ref)
for j in range(nr) : d+= ' '+fe(ref[j])+' '+fe(ewert[j])+'\n'
print(d); return d
def computedcase(subset, dbg=0, big=False) : # subsets: A B C D E F G
# case F,G needs w on x=-100..100,y=0..100, later node[-x,-y]=node[x,y]
def vscale(f,v) : # scale vector v
n=len(v); u=[0]*n
for i in range(n) : u[i]=f*v[i]
return u
nx=100; ny=100; w=[[]]*(nx+1)
nodiag=('EFG'.find(subset)>=0); xs= (-1) if(subset=='G') else 1
sign=(-1) if('BD'.find(subset)>=0) else 1 # EFG no x=y mirror!
rm,ixt,pars= ritzmatrix(subset,0,True,big)
# rm= zu diagonalisierende Operatormatrix
## eval,evec,rank= eigenrun(0,rm,dbg=1) # obsolete, use eigenloop
eval,evec= eigenloop(rm)
print(fdump(' eigenvals',eval))
ref,z = platewave(subset, 0,0, pars) # reference eigenval spectrum
if (subset!='G') and False : # deaktivierter Log
s= dumpmodes(subset,ixt, eval,evec,ref)
f=open('ritzmodes.txt','a'); f.write(s+'\n'); f.close()
for ix in range(nx+1) :
x= ix/(nx+0.0); w[ix]=[[]]*(ny+1); ymax= (nx+1) if nodiag else (ix+1)
for iy in range(ymax) :
y= iy/(ny+0.0)
w[ix][iy] = eigenwave(evec,ixt, xs*x,y,pars)
if not nodiag : # x=y mirror
for ix in range(nx+1) :
for iy in range(ix) : w[iy][ix]= vscale(sign,w[ix][iy]) # Symm-Achse x=y
print('Wellenwerte: '+str([len(w),len(w[0]),len(w[0][0])]))
if ('ABCDE'.find(subset)>=0) :
curves=knotenbild(subset, w, debug=dbg, version='2')
if (subset=='F')or(subset=='G') :
curves=knotenbild(subset, w, debug=3, version='2')
return curves
def ritzplots3() : # text Resultate, subset A B C D E F
f=open('ritzmodes.txt','w')
for subset in ['A','B','C','D','E','F'] :
rm,ixt,pars= ritzmatrix(subset,0,True,big=True)
ref,z = platewave(subset, 0,0, pars) # reference eigenval spectrum
eval,evec= eigenloop(rm); norm=eigentest(rm,eval,evec)
print('Subset='+subset+' Quality Eigenvals='+str(norm))
s= dumpmodes(subset,ixt, eval,evec,ref); f.write(s+'\n')
f.close()
def chladniquadrat() :
for c in ['A','B','C','D','E'] : simplecase(c, dbg=2)
specialcase() # diagonale 'F'-Serie (E bis)
from time import time
def ritzplots() : # Urversion
t0=time()
#simplecase('A',dbg=2)
chladniquadrat()
t1=time()
print('Laufzeit (sek)= '+str(int(t1-t0)))
def ritzplots2() : # Neuberechnung
t0=time()
computedcase('A', dbg=2, big=True)
computedcase('B', dbg=2, big=True)
computedcase('C', dbg=2, big=True)
computedcase('D', dbg=2, big=True)
computedcase('E', dbg=2, big=True)
specialcase(compute=2,version='2') # Serie 'F'
t1=time()
print('Laufzeit (sek)= '+str(int(t1-t0)))
def startdialog() :
try : s=input('Chladni-Figuren Urversion(1) Neurechnung(2) Eigenwerte(3)?')
except : s='???'
s= str(s)+'?'
if s[0]=='1' : ritzplots()
elif s[0]=='2' : ritzplots2()
elif s[0]=='3' : ritzplots3()
else : print('Eingabefehler.')
startdialog()
Listing chladni.py
Bearbeiten# -*- coding: utf-8 -*-
from math import sqrt
from ritz import eigenloop, eigentest, knotenbild, ritzmatrix
# def knotenbild(subset, w, debug=0, version='1') : # Figur aus Knotenlinien
# debug=2 macht Kurven mit Namen 'single'+subset+version
def fe(x) : return ('%11d' % x) if(type(x)==type(0)) else ('%+6.4e' % x)
def fdump(s,m) : # vector or matrix debug dump
n=len(m); ismat= (type(m[0])==type([])); t=''
if not ismat :
t='Vector '+s+' Dim='+str(n)+'\n'
for k in range(n) : t+=' '+fe(m[k])
else :
j=len(m[0]); t='Matrix '+s+' Dim='+str(n)+'*'+str(j)+'\n'
for i in range(n) :
for k in range(j) : t+=' '+fe(m[i][k])
if i<(n-1) : t+='\n'
return t
def norm(a) :
n=len(a); z=0.0
for i in range(n) : z += a[i]*a[i]
return sqrt(z)
def wsum(a,b,f,g) : # weighted sum of 2 vectors
n=len(a); c=[0]*n
for i in range(n) : c[i]= f*a[i]+g*b[i]
return c
def dpoly(p) : # Polynom-Ableitung
n=len(p); q=[0]*(n-1)
for i in range(n-1) : q[i]=p[i+1]*(i+1)
return q
def evpoly(p,x) : # Polynom-Auswertung
n=len(p); y=0; z=1
for i in range(n) : y+= z*p[i]; z *= x
return y
def tripleta(serie,max,dbg=0) : # serie=symmetry class. max= maximal order
# get Index triplets for a basis of symmetric linear combs
# i=k=0 excluded from function basis ?
m=1; grow=False; t=[] # m=target index + 1 with sign, t=triplet list.
if serie=='A' : i=1; k0=1; k=k0; sm=1 # sm= mirror symm
if serie=='B' : i=1; k0=1; k=k0; sm=-1
if serie=='C' : i=0; k0=0; k=k0; sm=1
if serie=='D' : i=0; k0=0; k=k0; sm=-1
if serie=='E' : i=0; k0=1; k=k0; sm=-1
if serie=='F' : i=0; k0=1; k=k0; sm=0
while (i<=max)and(k<=max) : # Kombis (0,0) (01) (10) sind immer steril ?
#if ((sm==1)and((i+k)>0)) or (i!=k) : t+= [[i,k,m]]; grow=True
## if ((i+k)>1) and ((sm==1) or (i!=k)) : t+= [[i,k,m]]; grow=True
if (sm>=0) or (i!=k) : t+= [[i,k,m]]; grow=True
if grow and (i!=k) and (sm!=0) : t+=[[k,i,sm*m]]
if grow : m+=1; grow=False
(i,k) = (i,k+2) if (k<i) else (i+2,k0)
if dbg==1 : print(str(t))
return t,m-1 # m is dim of target matrix to fill
### Funktional-Extrem mit orthonormalen Polynomen
def scalprod(p,q) : # Integral Polynome ueber -1...1
n=len(p); m=len(q); sp=0
for i in range(n) :
for k in range(m) :
j=i+k; z=(0 if(j%2>0) else 2.0/(j+1.0)); sp+= p[i]*q[k]*z
return sp
def spmatrix(p) : # debug
n=len(p); m=[[]]*n
for i in range(n) :
m[i]=[0]*n
for k in range(n) : m[i][k]= scalprod(p[i],p[k])
print(fdump('spmatrix',m))
def derivprod(p,q,np,nq) : # product of derivative (p,np) with (q,nq)
for i in range(np) : p= dpoly(p)
for i in range(nq): q= dpoly(q)
return scalprod(p,q) # needed derivprod(p,q,1,1), (p,q,2,0), (p,q,2,2)
def orthonormset(n,sprod,mode=0) :
# make polynomials p[0] to p[n], wrt scalar product sprod()
def normalize(p) :
z=sprod(p,p); q= 1.0/sqrt(z)
for j in range(len(p)) : p[j] *= q
n1=n+1; p=[[]]*n1; sp=[0]*n1
for i in range(n1) :
if (mode==0)or(i<2) : p[i]=[0]*(i+1); p[i][i]=1
elif mode==1 : # bigger polys with @boundary 2nd, 3rd deriv to vanish
ni=i+4; p[i]=[0]*(ni+1); p[i][i]=1; z=i*(i-1)
p[i][i+2]=-2*z/(i+1.0)/(i+2.0); p[i][i+4]= z/(i+3.0)/(i+4.0)
normalize(p[0])
for i in range(1,n1) : # orthogonalize p[i] wrt to subset {0...i-1}
for k in range(i) : sp[k]= sprod(p[i],p[k])
for k in range(i) :
for j in range(len(p[k])) : p[i][j]-= sp[k]*p[k][j]
normalize(p[i])
# for i in range(n1) : print(fdump('debug',p[i]))# are even/odd alternating
return p
def potentialmatrix(serie, max, nu, mode=0, dbg=0) :
# matrices on ortho poly sets
def auxmatrices(p) :
n=len(p); dp00=[[]]*n; dp11=[[]]*n; dp20=[[]]*n; dp22=[[]]*n
for i in range(n) :
dp00[i]=[0]*n; dp11[i]=[0]*n; dp20[i]=[0]*n; dp22[i]=[0]*n
for k in range(n) :
dp00[i][k]= scalprod(p[i],p[k])
dp11[i][k]= derivprod(p[i],p[k],1,1)
dp20[i][k]= derivprod(p[i],p[k],2,0)
dp22[i][k]= derivprod(p[i],p[k],2,2)
return dp00, dp11, dp20, dp22
#
def vorthomatrix(n,m,p,q,nu,dp00,dp11,dp20,dp22) : # pot. on orthog (nm|pq)
# V(f,g)=fxx*gxx+fyy*gyy+nu(fxx*gyy+fyy*gxx)+2(1-nu)(fxy*gxy)
a= dp22[n][p]*dp00[m][q] # derivpr(p[n],p[p],2,2)
b= dp22[m][q]*dp00[n][p] # derivpr(p[m],p[q],2,2)
c= dp11[n][p]*dp11[m][q] # derivpr(p[n],p[p],1,1)*derivpr(p[m],p[q],1,1)
d= dp20[n][p]*dp20[q][m] # derivpr(p[n],p[p],2,0)*derivpr(p[m],p[q],0,2)
e= dp20[p][n]*dp20[m][q] # derivpr(p[n],p[p],0,2)*derivpr(p[m],p[q],2,0)
f=1;g=2 #dbg
v= (a+b) + f*nu*(d+e) + g*(1-nu)*c; s=dp00[n][p]*dp00[m][q] #s=scalprod
return v,s
#
pol= orthonormset(max,scalprod,mode=mode)
dp00, dp11, dp20, dp22 = auxmatrices(pol)
ixtrip,mmat= tripleta(serie,max,dbg=dbg); nix=len(ixtrip)
pm=[[]]*mmat; nm=[[]]*mmat; sm=[[]]*mmat;
for h in range(mmat) : pm[h]=[0]*mmat; nm[h]=[0]*mmat; sm[h]=[0]*mmat
for i in range(nix) : # add to row xi
mi,ni,xi= tuple(ixtrip[i]); (si,xi) =(1,xi-1) if (xi>0) else (-1,-xi-1)
for j in range(nix) : # addto column xj
mj,nj,xj= tuple(ixtrip[j]); (sj,xj) =(1,xj-1) if (xj>0) else (-1,-xj-1)
v,s= vorthomatrix(mi,ni, mj,nj, nu, dp00,dp11,dp20,dp22)
pm[xi][xj]+= si*sj*v; sm[xi][xj]+= si*sj*s; nm[xi][xj]+=1
for i in range(mmat) : # normalize
for j in range(mmat) :
if nm[i][j]>1 : z=1.0/sqrt(nm[i][j]); pm[i][j]*=z; sm[i][j]*=z
if dbg==1 : print(fdump('pm',pm)) # potential energy matrix, bug serie E?
if dbg==1 : print(fdump('sm',sm)) # scalar prod matrix seems correct.
return pm, ixtrip, pol
def showsorted(ser,v,j) : # lowest positive values in vector v
n=len(v); ix=[-1]*n; iy=[0]*n; ival=0; s='Serie '+ser+':'
bound=1e-12; vsort=[0]*j
for i in range(n) :
vmin=1e33; kmin=-1
for k in range(n) :
if (ix[k]<0) and (v[k]<vmin) and (v[k]>bound) : vmin=v[k]; kmin=k
iy[i]=kmin; ix[kmin]=ival; ival+=1
for i in range(min(j,n)) : s+=' '+fe(v[iy[i]]); vsort[i]=v[iy[i]]
return vsort,s
def wellenbilder(subset,eval,evec,ixt,pol) : # Eigenkombinationen aus Polynomen
# indizes der 5 kleinsten positiven Eigenwerte
def auswahl(eva,nmax) :
m=len(eva); nmax=min(nmax,m); ix=[-1]*nmax; use=[0]*m; bound=1e-12
for i in range(nmax) :
imin=-1; emin=1e33
for k in range(m) :
if (use[k]==0)and(eva[k]>bound)and(eva[k]<emin): (imin,emin)=(k,eva[k])
ix[i]=imin; bound=eva[imin]; use[imin]=1
print('Auswahl:'+str(ix))
return ix
#
def eigenwave(iew, ev, ixt,x,y,pol) : # Wellen aus Vektoren ev in Basis ixt
np=len(pol); u=[0]*np; v=[0]*np; new=len(iew)
for i in range(np) : u[i]= evpoly(pol[i],x); v[i]= evpoly(pol[i],y)
bw=[0.0]*50; kmax=0; n=len(ixt) # bw=basiswellen, w=eigenwellen
for i in range(n) : # Basis-Wellenwerte
j,k,nd= tuple(ixt[i]); (sgn,ix)= (1,nd-1) if (nd>0) else (-1,-nd-1)
bw[ix]+= sgn*u[j]*v[k]; kmax=max(kmax,ix)
m=kmax+1; w=[0.0]*new # linearkombiniert zu Eigenwellen
for h in range(new) :
i=iew[h]
for j in range(m) : w[h]+= ev[i][j]*bw[j]
return w # m Ergebnisse am Punkt (x,y)
#
def vscale(f,v) : # scale vector v
n=len(v); u=[0]*n
for i in range(n) : u[i]=f*v[i]
return u
#
iew= auswahl(eval,5); new=len(iew)
nx=100; ny=100; w=[[]]*(nx+1); xs=1
for ix in range(nx+1) :
x= ix/(nx+0.0); w[ix]=[[]]*(ny+1)
for iy in range(ny+1) :
y= iy/(ny+0.0); w[ix][iy] = eigenwave(iew, evec,ixt, xs*x,y,pol)
print('Wellenwerte: '+str([len(w),len(w[0]),len(w[0][0])]))
knotenbild(subset, w, debug=2, version='3') # ritz.knotenbild
return w
def tryorthopolyn(maxi, version, dbg=0) :
print('Orthogonale Polynome maxi='+str(maxi)+' version='+str(version))
nu=0.225;
for serie in ['A','B','C','D'] : ## ,'E','F'] :
ma,ixt,pol= potentialmatrix(serie,maxi,nu, mode=version, dbg=dbg)
eva,xa= eigenloop(ma)
eigentest(ma,eva,xa, dbg=dbg)
vsort, ssort= showsorted(serie,eva,5); print(ssort)
wellenbilder(serie,eva,xa,ixt,pol) # knotenbilder
def spektrumvergleich(maxi) :
nu=0.225; version=0; print('Spektrum Balkenbasis Polynombasis')
for serie in ['A','B','C','D'] :
ma,ixt,pol= potentialmatrix(serie,maxi,nu, mode=version, dbg=0)
eval,xa= eigenloop(ma)
vsort,ssort = showsorted(serie,eval,5)
rm,ixr,pars= ritzmatrix(serie,0,True,big=True)
reval,revec= eigenloop(rm); # print(str(reval[:5]))
for k in range(5) :
s=serie+str(k+1)+' '+fe(reval[k])+' '+fe(vsort[k]); print(s)
def startdialog() :
s=str(input('Polynombasis: Chladni-Figuren(1) Eigenwert-Vergleich(2)?'))+'?'
if s[0]=='1' : tryorthopolyn(15,0) # production run
elif s[0]=='2' : spektrumvergleich(15)
else : print('Eingabefehler.')
startdialog()