Zum Inhalt springen

Till Eulenspiegels lustige Serie/ Schwingende Platten/ Skripte

Aus Wikibooks

Kurzbeschreibung

[Bearbeiten]

Die 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()