p:=matrix([x,y,z]):pt:=linalg::transpose(p)
Quadriken 3d Ellipsoid
Prof. Dr. Dörte Haftendorn: Mathematik mit MuPAD 4, Juni 07 Update 5.07.07
Web: https://mathe.web.leuphana.de www.mathematik-verstehen.de ######################################################
Rückwärts aufgestellt, siehe Datei Konstruktion-rückwärts
p:=matrix([x,y,z]):pt:=linalg::transpose(p)
quadrik:=matrix([[7*x^2 - 4*x*y - 2*x*z - 70*x + 7*y^2 - 2*y*z + 38*y + 4*z^2 + 28*z + 202]])
quadrikp:=plot::Implicit3d(quadrik[1],x=-2..7,y=-5..0,z=-7..0,FillColor=[0,1,0,0.8], Scaling=Constrained):
plot(quadrikp)
A und a passend aufstellen
A:=matrix([[7,-2,-1],[-2,7,-1],[-1,-1,4]]);
a:=matrix([-70,+38,28]): at:=linalg::transpose(a);
d:=202;
expand(pt*A*p+at*p+d); //Probe, ob man A,a und d richtig hat
%-quadrik
hier muss 0 herauskommen
Hauptachtsentransformation
E3:=matrix([[1,0,0],[0,1,0],[0,0,1]])
evli:=linalg::eigenvectors(A) //Probe, was MuPAD liefert
Eigenwerte und Eigenvektoren
ew1 :=evli[1][1]; ew2 :=evli[2][1]; ew3 :=evli[3][1];
ev1:=evli[1][3][1];
ev2:=evli[2][3][1];
ev3:=evli[3][3][1];
ev1n:=linalg::normalize(ev1):
ev2n:=linalg::normalize(ev2):
ev3n:=linalg::normalize(ev3):
P:=ev1n.ev2n.ev3n: Pt:=linalg::transpose(P);
Vektorschreibweise für die Abbildung und die
Quadrikgleichungen, die sich durch Einsetzen ergeben:
bzw. in 3D
quastrich:=Simplify(pt*Pt*A*P*p+at*P*p+d)
quastrichp:=plot::Implicit3d(quastrich[1],x=-5..0,y=-5..4, z=-7..-2,FillColor=[1,1,0,1]):
plot(quastrichp /*,Op, tp */, Scaling=Constrained):
quastrich
Die Arbeitsweise ist dieselbe wie bei der Herstellung der Scheitelform einer Parabel.
Hier durch Hinsehen:
xterm:=hold(6*(x+5/3*sqrt(3))^2);expand(xterm);
yterm:=hold(3*(y+2/3*sqrt(6))^2);expand(yterm);
zterm:=hold(9*(z+3*sqrt(2))^2);expand(zterm);
Also
quastrichK:=xterm+yterm+zterm-18 //-220+202
quastrich-expand(quastrichK)
hier muss 0 herauskommen
Letzter Teil der Hauptachsentransformation ist die Translation t
t:=matrix([5/3*sqrt(3), 2/3*sqrt(6),3*sqrt(2)]);
tt:=linalg::transpose(t):
also
Das ergibt:
quaH:=Simplify(expand((pt-tt)*Pt*A*P*(p-t)+at*P*(p-t)))+202;
quadrikH:=ew1*x^2+ew2*y^2+ew3*z^2-18
Angabe der Gleichung in der üblichen Form:
hold(x^2/3+y^2/6+z^2/2=1)
quadrikHp:=plot::Implicit3d(quadrikH,x=-2..2,y=-3..3, z=-2..2,
Scaling=Constrained):
plot(%)
Bestimmung des ursprünglichen Mittelpunktes:
nun in 3D
m:=Simplify(P*(-t))
Dieses ist der alte Mittelpunkt.
Bestimmung des Urbildes des rechten Hauptscheitels:
Mit ev1 rechts ist hier der normierte 1. Eigenvektor gemeint. Alles jetzt in 3D
r:=sqrt(3): // Halbachse in (späterer) x-Richtung
rur:=r*ev1n+m;
rurstrich:=Pt*rur;
mp:=plot::Point3d(m,PointSize=2, PointColor=[0,1,0]):
Op:=plot::Point3d([0,0,0], PointSize=2, PointColor=[1,0,0]):
rurp:=plot::Point3d(rur, PointSize=2, PointColor=[0,1,1]):
rursp:=plot::Point3d(rurstrich, PointSize=2, PointColor=[1,0.5,0]):
rp:=plot::Point3d([r,0,0],PointSize=2, PointColor=[0,0,1]):
ev1urp:=plot::Arrow3d(m,m+3*ev1):
ev2urp:=plot::Arrow3d(m,m+3*ev2):
ev3urp:=plot::Arrow3d(m,m+3*ev3):
tp:=plot::Arrow3d(-t,[0,0,0],LineColor=[1,0,0]):
plot(quadrikp,mp,Op,rurp,tp,
ev1urp,ev2urp,ev3urp,PointSize=2,Scaling=Constrained);
Test ob es eine Drehung allein oder gefolgt von Spiegelung ist. det=1 reine Drehung
linalg::det(Pt)
plot(quadrikp,quastrichp, quadrikHp,mp,Op,tp,rp,rurp,rursp,
ev1urp,ev2urp,ev3urp,Axes=Origin,Scaling=Constrained)
####################################################
Weitere Untersuchungen,
die Abbildung wird von Pt geleistet
Pt
Suche nach der Drehachse, sie muss die Eigenrichtung der Drehmatrix sein.
evliP:=linalg::eigenvectors(Pt)
float(evliP)
Wie erwartet ist nur ein Eigenwert reell und zwar ist er 1.
Drehachsenrichtung:
dreh:=evliP[1][3][1]
Die Idee, dass die Normale auf der Ebene durch 0, M und M' die Drehachse sei, ist falsch:
mtt:=linalg::crossProduct(-t,m); float(mtt)
Obiger Vektor stimmt nicht mit den Normalenvektor überein.
drehp:=plot::Arrow3d(5*dreh, LineColor=[1,0,1])
Einzeichnen der Drehachse:
plot(quadrikp,quastrichp, quadrikHp,mp,Op,tp,rp,rurp,rursp,
ev1urp,ev2urp,ev3urp,drehp,Axes=Origin,Scaling=Constrained)
Die besondere Herausfordung, diese Drehung animiert zu gestalten, habe ich noch
nicht umgesetzt.
Zunächst muss ich herausbekommen, in welchen Ebenen sich die Punkte jeweils bewegen.
Hessesche Normalenform:
float(linalg::scalarProduct( (p-m),dreh)=0);
float(linalg::scalarProduct( (p+t),dreh)=0)
Hier sieht man, dass tatsächlich m und mstrich=-t in dieser Ebene liegen.
Wo durchstößt die Drehachse diese Ebene?
Dieser Durchstoßpunkt muss das k-fache des Drehachsenvektors sein. )
1.303225373^2*k - 0.1109881895^2*k + 1.0^2*k + 6.050403504 = 0.0;
solve(%)
delete k:
Der Durchstoßpunkt ist also
k:= -2.252504773: dm:=float(k*dreh)
Bestimmung des Drehwinkels
m-dm, -t-dm, linalg::scalarProduct(m-dm,-t-dm)
cosi:=linalg::scalarProduct(m-dm,-t-dm)/ sqrt(linalg::scalarProduct(m-dm,m-dm)*
linalg::scalarProduct(-t-dm,-t-dm))
float(arccos(cosi))
################################################################
Das mache ich jetzt zur Probe auch für den rechten Hauptscheitel
rurs:=rurstrich:
float(linalg::scalarProduct( (p-rur),dreh)=0);
float(linalg::scalarProduct( (p-rurs),dreh)=0)
Hier sieht man, dass tatsächlich rur und rurstrich=rurs in derselben Ebene liegen.
Wo durchstößt die Drehachse diese Ebene?
Dieser Durchstoßpunkt muss das k-fache des Drehachsenvektors sein. )
delete kk:
1.303225373^2*kk - 0.1109881895^2*kk + 1.0^2*kk +6.242640687 = 0.0;
solve(%)
Der Durchstoßpunkt ist also
kk:= -2.324072756: rm:=float(kk*dreh)
rur-rm, rurs-rm, linalg::scalarProduct(rur-rm,rurs-rm)
cosi:=linalg::scalarProduct(rur-rm,rurs-rm)/
sqrt(linalg::scalarProduct(rur-rm,rur-rm)*
linalg::scalarProduct(rurs-rm,rurs-rm))
wradm:=float(arccos(cosi));
wgrdm:=float(wradm/PI*180);
Das sind dieselben Werte wie für m und mstrich
wradm;
wgrdm;
Damit handelt es sich also um eine Drehung um 125,75.. Grad.