Auswahl Aktion im Forum Suche Kontakt Für Mitglieder Mathematisch für Anfänger Wer ist Online | |
| Autor |
 Projektive Ebene mit maple16 Projekt |
|
calabi-yau
Senior 
Dabei seit: 15.04.2006
Mitteilungen: 901
Aus:
 |      Themenstart: 2012-08-22 23:27
|
Hallo,
ich möchte die Projektive Ebene ähnlich wie hier mit Papierstreifen nachbauen. Mathematisch habe ich die Sache schon durchgedacht, aber ich möchte keinesfalls das alles per hand durchrechnen und auch das plotten wird ohne Computer schwer. Ich möchte daher maple(16) verwenden, allerdings reichen meine rundimentären Kenntnisse nicht mal für den ersten Schritt aus und ich möchte auch nicht erst ein Buch durchwälzen. Ich möchte es daher mal mit dem Forum probieren und beginne mit dem ersten Schritt:
Zuerst soll folgende Funktion
und dann die Funktion
in Maple bekanntgemacht werden und zur Kontrolle soll das Bild von B geplottet werden. Kann mir jemand da eine Starthilfe geben?
|
 Profil
 Quote
 Link |
Ollie
Senior 
Dabei seit: 03.05.2003
Mitteilungen: 5866
Aus: Aachen
 | Beitrag No.1, eingetragen 2012-08-23 01:32
|
Hallo, hier ein Ansatz.
Maple restart;with(plots):
1/(z^3-z^(-3)+sqrt(5))*(<I*(z^2-z^(-2)),z^2+z^(-2),2*I/3*(z^3+z^(-3))>);
subs(z=x+I*y,%);
Re(%)assuming real:M:=%+<0,0,1/2>;
M0:=LinearAlgebra:-VectorNorm(M,2):
B:=M/M0^2:
B1:=subs(y=1,B); # Schnittebene y=1 als Beispiel, k.A. ob das sinnvoll ist
spacecurve([B1(1),B1(2),B1(3)],x=-1e5..1e5,colour=black,axes=normal); |
Gruß
[ Nachricht wurde editiert von Ollie am 23.08.2012 01:41:55 ]
|
Profil
Quote
Link |
calabi-yau
Senior
Dabei seit: 15.04.2006
Mitteilungen: 901
Aus:
| Beitrag No.2, vom Themenstarter, eingetragen 2012-08-23 14:54
|
Hi Ollie, vielen Dank für deine Hilfe. Das sieht schon mal gut aus. ;)
Also folgendes steht jetzt bei mir im Worksheet:
Maple restart;
with(plots):with(LinearAlgebra):
M:= Re( 1/(z^3-z^(-3)+sqrt(5)) * <I*(z^2-z^(-2)),z^2+z^(-2),2*I/3*(z^3+z^(-3))> ) + <0,0,2/5>:
B:=M/VectorNorm(M,2)^2:
BReal:=subs(z=x+I*y,B):
#plot3d([BReal(1),BReal(2),BReal(3)],x=-1..1,y=-1..1,grid=[80, 80]);
c:=I*t:
Bc:=subs(z=c,B):
#spacecurve([Bc(1),Bc(2),Bc(3)],t=0..1,color=black,axes=normal);
BcL:=int(VectorNorm(Bc,2),t=0..s):
#plot(BcL,s=0..1); |
Also ich hab die Kurve Bc (das wird eine von vielen) und was ich nun machen will ist diese auf Bogenlänge zu parametrisieren. Das erste was mir auffällt ist, das bereits jetzt dies schon etwas rechenintensiv ist, dabei würde es mir numerisch absolut genügen. Kann man das durch einen Befehl festlegen oder muss man da ganz anders vorgehen?
Ansonsten habe ich jetzt die Längenfunktion BcL und was ich jetzt machen muss, ist diese Expression BcL nach s aufzulösen, also einfach die Funktion zu invertieren. Wie gehe ich da vor?
Edit: Ich merke gerade, das Integral symbolisch auszurechnen dauert nicht lange, nur der plot, also fällt diese Frage schon mal weg. ;)
[ Nachricht wurde editiert von calabi-yau am 24.08.2012 12:25:02 ]
|
Profil
Quote
Link |
Ollie
Senior
Dabei seit: 03.05.2003
Mitteilungen: 5866
Aus: Aachen
| Beitrag No.3, eingetragen 2012-08-24 13:16
|
zum Plot: Ersetzte int durch Int, dann geht das Plotten in Sekunden.
Am Schluß habe ich noch eine Zeile angehängt, die numerisch die obere Integrationsgrenze bei gegebener Bogenlänge berechnet.
Ach ja: Durch Einfügen von assume(t::real); wurde BcL übersichtlicher.
Gruß
Maple restart; assume(t::real);
with(plots):with(LinearAlgebra):
M:= Re( 1/(z^3-z^(-3)+sqrt(5)) * <I*(z^2-z^(-2)),z^2+z^(-2),2*I/3*(z^3+z^(-3))> ) + <0,0,2/5>:
B:=M/VectorNorm(M,2)^2:
BReal:=subs(z=x+I*y,B):
#plot3d([BReal(1),BReal(2),BReal(3)],x=-1..1,y=-1..1,grid=[80, 80]);
c:=I*t:
Bc:=subs(z=c,B):
#spacecurve([Bc(1),Bc(2),Bc(3)],t=0..1,color=black,axes=normal);
BcL:=Int(VectorNorm(Bc,2),t=0..s);
plot(BcL,s=0..1);
BcLn:=proc(L) fsolve(L=Int(VectorNorm(Bc,2),t=0..s),s); end proc;
BcLn(0.5); |
|
Profil
Quote
Link |
calabi-yau
Senior
Dabei seit: 15.04.2006
Mitteilungen: 901
Aus:
| Beitrag No.4, vom Themenstarter, eingetragen 2012-08-24 15:09
|
Hi, vielen Dank, wir nähern uns in kleinen Schritten einem guten Anfang. ;)
Nun sieht der Code bei mir so aus, nur funktionieren die letzten beiden Codezeilen nicht so wie ich es mir vorstelle. Den plot will ich mal wieder zur Kontrolle machen, das Ableiten sollte das normierte Tangentialvektorfeld zu Kurve geben.
Maple restart; assume(t::real);
with(plots):with(LinearAlgebra):
M:= Re( 1/(z^3-z^(-3)+sqrt(5)) * <I*(z^2-z^(-2)),z^2+z^(-2),2*I/3*(z^3+z^(-3))> ) + <0,0,2/5>:
B:=M/VectorNorm(M,2)^2:
BReal:=subs(z=x+I*y,B):
#plot3d([BReal(1),BReal(2),BReal(3)],x=-1..1,y=-1..1,grid=[80, 80]);
c:=I*t:
Bc:=subs(z=c,B):
#spacecurve([Bc(1),Bc(2),Bc(3)],t=0..1,color=black,axes=normal);
BcL:=Int(VectorNorm(Bc,2),t=0..s):
l1:=subs(s=0.2,BcL):
l2:=subs(s=0.4,BcL):
l3:=subs(s=0.6,BcL):
l4:=subs(s=0.8,BcL):
lmax:=subs(s=1,BcL):
#plot(BcL,s=0..1);
BcLi:=proc(L) fsolve(L=BcL,s); end proc:
cc:=unapply(Bc,t)@BcLi:
spacecurve(cc(s),s=0..lmax); #?!
ccD:=D(cc); #?! |
|
Profil
Quote
Link |
Ollie
Senior
Dabei seit: 03.05.2003
Mitteilungen: 5866
Aus: Aachen
| Beitrag No.5, eingetragen 2012-08-24 22:22
|
Ich habe die Prozedur cc jetzt nicht direkt an spacecurve übergeben, sondern diese auf eine diskrete Punkteliste gemapt und dann letztere gespacecurved. Was soll eigentlich das Endergebnis unserer Bemühungen sein? :-)
Maple restart; assume(t::real,c::real);
with(plots):with(LinearAlgebra):
M:= Re( 1/(z^3-z^(-3)+sqrt(5)) * <I*(z^2-z^(-2)),z^2+z^(-2),2*I/3*(z^3+z^(-3))> ) + <0,0,2/5>:
B:=normal(M/VectorNorm(M,2)^2):
BReal:=subs(z=x+I*y,B):
#plot3d([BReal(1),BReal(2),BReal(3)],x=-1..1,y=-1..1,grid=[80, 80]);
c:=I*t:
Bc:=subs(z=c,B):Bc:=simplify(Bc,symbolic) assuming real;
#spacecurve([Bc(1),Bc(2),Bc(3)],t=0..1,color=black,axes=normal);
BcL:=Int(VectorNorm(Bc,2),t=0..s):
l1:=subs(s=0.2,BcL):
l2:=subs(s=0.4,BcL):
l3:=subs(s=0.6,BcL):
l4:=subs(s=0.8,BcL):
lmax:=subs(s=1,BcL):
#plot(BcL,s=0..1);
BcLi:=proc(L) fsolve(L=BcL,s); end proc:
cc1:=unapply(Bc,t):
cc:=cc1@BcLi:
seq(n/10,n=-20..20):
map(cc,[%]):
S:=map(convert,%,list):
spacecurve(S,colour=black,axes=normal,style=line);
|
|
Profil
Quote
Link |
calabi-yau
Senior
Dabei seit: 15.04.2006
Mitteilungen: 901
Aus:
| Beitrag No.6, vom Themenstarter, eingetragen 2012-08-24 23:26
|
Das Endergebnis ist ein Nachbau der reellen projectiven Ebene als Immersion in den R^3 nach Werner Boy (Die Funktion B oben stammt genauer von Bryant/Kusner) wie sie in Oberwolfach steht (siehe Link oben). Dazu nimmt man einige Parameterlinien und plottet sie auf Papier. Genauer nimmt man eine Parameterlinie c wie oben. Nun muss man eine Kurve auf dem flachen Blatt Papier so plotten, dass deren (signierte) Krümmung im Flachen gleich der intrinsischen (signierten) Krümmung von c auf der Boyschen Fläche ist. Dies hat den Effekt, dass wenn man die geplottete Linie zu einem Streifen aufdickt, dieser dann tangential zur Fläche verläuft, wenn man alle Streifen zusammenbaut. Dazu helfen Punktierungen auf den Streifen, die die Schnittpunkte der Streifen markieren. Wie man das im genauen macht, sehen wir, wenn wir im Maple Projekt weitermachen. ;)
Hast du eine Lösung für die Ableitung von cc? Ich meine ich sehe das Problem: Bereits BcLi ist nur numerisch gegeben, kann maple keine numerischen Funktionen ableiten?
PS: Die reelle projective Ebene ist die eindeutige geschlossene, nichtorientierte, irreducible (bzgl. zusammenhängender Summe #) 2-Mannigfaltigkeit und neben der Sphäre und dem Torus die dritte Fläche, mit der sich alle geschlossenen 2-Mannigfaltigkeiten klassifizieren lassen. Die Immersion oben minimiert außerdem das Willmore-Funktional (minimale Elastizitätsenergie). Ein wunderbares Model um Laien etwas Topologie näherzubringen.
|
Profil
Quote
Link |
Ollie
Senior
Dabei seit: 03.05.2003
Mitteilungen: 5866
Aus: Aachen
| Beitrag No.7, eingetragen 2012-08-24 23:57
|
Aha.Topologisch unbeleckt wie ich bin, nehme ich das erstmal zur Kenntnis und zeige noch ein Beispiel zur numerischen Differentation in maple.
Maple f := proc (t) t^2 end proc;
fdiff(f, [1], [3]); |
Die erste eckige Klammer gibt die Differentationsvariable(n) an, die zweite den Auswertepunkt.
|
Profil
Quote
Link |
calabi-yau
Senior
Dabei seit: 15.04.2006
Mitteilungen: 901
Aus:
| Beitrag No.8, vom Themenstarter, eingetragen 2012-08-25 01:46
|
Edit: siehe unten
[ Nachricht wurde editiert von calabi-yau am 26.08.2012 01:53:23 ]
|
Profil
Quote
Link |
calabi-yau
Senior 
Dabei seit: 15.04.2006
Mitteilungen: 901
Aus:
| Beitrag No.9, vom Themenstarter, eingetragen 2012-08-25 19:03
|
Edit: siehe unten
[ Nachricht wurde editiert von calabi-yau am 26.08.2012 01:53:38 ]
|
Profil
Quote
Link |
calabi-yau
Senior
Dabei seit: 15.04.2006
Mitteilungen: 901
Aus:
| Beitrag No.10, vom Themenstarter, eingetragen 2012-08-26 01:52
|
also ich habe nun eine möglichkeit gefunden, die parametrisierung nach bogenlänge zu vermeiden. zur besseren übersicht habe ich die 2 vorhergehenden (jetzt redundanten) posts gelöscht.
Maple restart; assume(t::real,x::real,y::real);
with(plots):with(LinearAlgebra):
M:= Re( 1/(z^3-z^(-3)+sqrt(5)) * <I*(z^2-z^(-2)),z^2+z^(-2),2*I/3*(z^3+z^(-3))> ) + <0,0,2/5>:
B:=M/VectorNorm(M,2)^2:
BReal:=subs(z=x+I*y,B): BReal:=simplify(BReal,symbolic) assuming real:
#plot3d([BReal(1),BReal(2),BReal(3)],x=-1..1,y=-1..1,grid=[80, 80]);
c:=I*t:
cReal:=<Re(c),Im(c)>:
k:=subs(x=cReal(1),y=cReal(2),BReal):
Dk:=<diff(k(1),t),diff(k(2),t),diff(k(3),t)>:
Dkn:=Dk/VectorNorm(Dk,2):
DxB:=<diff(BReal(1),x),diff(BReal(2),x),diff(BReal(3),x)>:
DyB:=<diff(BReal(1),y),diff(BReal(2),y),diff(BReal(3),y)>:
DxBc:=subs(x=cReal(1),y=cReal(2),DxB):
DyBc:=subs(x=cReal(1),y=cReal(2),DyB):
N:=CrossProduct(DxBc,DyBc):
Nn:=N/VectorNorm(N,2):
q:=Matrix([[1+Nn(1)*I,Nn(2)+Nn(3)*I],[-Nn(2)+Nn(3)*I,1-Nn(1)*I]]):
x:=Matrix([[Dkn(1)*I,Dkn(2)+Dkn(3)*I],[-Dkn(2)+Dkn(3)*I,-Dkn(1)*I]]):
vv:=MatrixMatrixMultiply(q,MatrixMatrixMultiply(x,MatrixInverse(q))):
v:=<Im(vv[1,1]),Re(vv[1,2]),Im(vv[1,2])>:
DDkn:=<diff(Dkn(1),t),diff(Dkn(2),t),diff(Dkn(3),t)>:
w:=DDkn/VectorNorm(Dk,2):
f:=v(1)*w(1)+v(2)*w(2)+v(3)*w(3):
#evalf(subs(t=0.1,f));
#plot(f,t=0..1,numpoints=10);
chi:=Int(f*VectorNorm(Dk,2),t=0..a):
#evalf(subs(a=0.5,chi));
res1:=Int(cos(chi)*subs(t=a,VectorNorm(Dk,2)),a=0..s):
res2:=Int(sin(chi)*subs(t=a,VectorNorm(Dk,2)),a=0..s):
evalf(subs(s=0.5,res1)); |
Das wäre das Endergebnis (für einen Streifen), allerdings schafft mein PC es nicht, den letzten Befehl in annehmbarer Zeit zu berechnen. Was könnte man vereinfachen?
|
Profil
Quote
Link |
Ollie
Senior
Dabei seit: 03.05.2003
Mitteilungen: 5866
Aus: Aachen
| Beitrag No.11, eingetragen 2012-08-26 11:10
|
Kein Wunder, der Teilintegrand f ist Dutzende worksheet-Seiten lang.
Nach symbolischen Berechnungen sollte man immer die Teilausdrücke untersuchen und ggfs. vereinfachen lassen. Keep it simple.
Die letzte Zeile gibt neben dem Ergebnis auch die Berechnungszeit des subs-Befehls aus.
Maple restart; assume(t::real,x::real,y::real);
with(plots):with(LinearAlgebra):
M:= Re( 1/(z^3-z^(-3)+sqrt(5)) * <I*(z^2-z^(-2)),z^2+z^(-2),2*I/3*(z^3+z^(-3))> ) + <0,0,2/5>:
B:=M/VectorNorm(M,2)^2:
BReal:=subs(z=x+I*y,B): BReal:=simplify(BReal,symbolic) assuming real:
#plot3d([BReal(1),BReal(2),BReal(3)],x=-1..1,y=-1..1,grid=[80, 80]);
c:=I*t:
cReal:=<Re(c),Im(c)>:
k:=subs(x=cReal(1),y=cReal(2),BReal):
Dk:=<diff(k(1),t),diff(k(2),t),diff(k(3),t)>:
Dkn:=Dk/VectorNorm(Dk,2):
DxB:=<diff(BReal(1),x),diff(BReal(2),x),diff(BReal(3),x)>:
DyB:=<diff(BReal(1),y),diff(BReal(2),y),diff(BReal(3),y)>:
DxBc:=subs(x=cReal(1),y=cReal(2),DxB):
DyBc:=subs(x=cReal(1),y=cReal(2),DyB):
N:=CrossProduct(DxBc,DyBc):
Nn:=N/VectorNorm(N,2):
q:=Matrix([[1+Nn(1)*I,Nn(2)+Nn(3)*I],[-Nn(2)+Nn(3)*I,1-Nn(1)*I]]):
x:=Matrix([[Dkn(1)*I,Dkn(2)+Dkn(3)*I],[-Dkn(2)+Dkn(3)*I,-Dkn(1)*I]]):
vv:=MatrixMatrixMultiply(q,MatrixMatrixMultiply(x,MatrixInverse(q))):
v:=<Im(vv[1,1]),Re(vv[1,2]),Im(vv[1,2])>:
DDkn:=<diff(Dkn(1),t),diff(Dkn(2),t),diff(Dkn(3),t)>:
w:=DDkn/VectorNorm(Dk,2):
f:=v(1)*w(1)+v(2)*w(2)+v(3)*w(3):
#evalf(subs(t=0.1,f));
#plot(f,t=0..1,numpoints=10);
f_simple:=simplify(f*VectorNorm(Dk,2)):
chi:=Int(f_simple,t=0..a):
#evalf(subs(a=0.5,chi));
VN:=simplify(VectorNorm(Dk,2)):
res1:=Int(cos(chi)*subs(t=a,VN),a=0..s):
#res2:=Int(sin(chi)*subs(t=a,VectorNorm(Dk,2)),a=0..s):
st:=time():evalf(subs(s=0.5,res1));time()-st;
|
|
Profil
Quote
Link |
calabi-yau
Senior
Dabei seit: 15.04.2006
Mitteilungen: 901
Aus:
| Beitrag No.12, vom Themenstarter, eingetragen 2012-08-26 15:11
|
Hi Ollie,
ja danke, das funktioniert. Zum Geschwindigkeitsvergleich wollte ich noch eine andere Methode ausprobieren, bei der ich die Umkehrfunktion der Längenfunktion doch wieder benutze, aber die fsolve Lösung ist irgendwie inkompatibel zum Rest:
Maple kL:=Int(DkNorm,t=0..a):
kLi:=proc(l) fsolve(l=kL,a); end proc:
fL:=proc(L) subs(t=kLi(L),f); end proc:
sys:=diff(theta(s),s)=fL(s), diff(b1(s),s)=cos(theta(s)), diff(b2(s),s)=sin(theta(s));
p:=dsolve({sys,theta(0)=0,b1(0)=0,b2(0)=0},{theta(s),b1(s),b2(s)},type=numeric):
test:=odeplot(p,[b1(s),b2(s)],0..5,numpoints=10):
display(test); |
Bei sys:=... kommt dann logischerweise
"Error, (in fsolve) s is in the equation, and is not solved for"
Wie packe ich die Umkehrfunktion denn richtig in die ODE?
PS: Ich fahre morgen in den Urlaub, evtl. kann ich dann nicht antworten, aber ich komme wieder. ;)
[ Nachricht wurde editiert von calabi-yau am 26.08.2012 15:11:32 ]
|
Profil
Quote
Link |
Ollie
Senior
Dabei seit: 03.05.2003
Mitteilungen: 5866
Aus: Aachen
| Beitrag No.13, eingetragen 2012-08-27 14:20
|
Hi, ich nehme an, dies ist nur ein Auschnitt aus einem anderen code als der bisherige.
Schau mal bei dsolve,numeric nach, unter dem Beispiel zur Option "known".
Gruß
|
Profil
Quote
Link |
calabi-yau
Senior
Dabei seit: 15.04.2006
Mitteilungen: 901
Aus:
| Beitrag No.14, vom Themenstarter, eingetragen 2012-09-12 21:47
|
Hallo Ollie,
sorry das ich solange nicht weitergeschrieben habe. Erst war ich im Urlaub und dann musste ich umziehen. Der Code sieht momentan folgendermaßen aus und funktioniert ganz gut:
Maple restart; assume(x::real,y::real,t::real);
with(plots):with(LinearAlgebra):
M:=Re(1/(z^3-z^(-3)+sqrt(5))*<I*(z^2-z^(-2)),z^2+z^(-2),2*I/3*(z^3+z^(-3))>)+<0,0,2/5>:
B:=M/VectorNorm(M,2)^2:
BReal:=subs(z=x+I*y,B): BReal:=simplify(BReal) assuming real:
#plot3d([BReal(1),BReal(2),BReal(3)],x=-1..1,y=-1..1,grid=[80, 80]);
c:=exp(I*(Pi/2))*t:
cReal:=<Re(c),Im(c)>:
k:=subs(x=cReal(1),y=cReal(2),BReal): k:=simplify(k) assuming real:
#spacecurve([k(1),k(2),k(3)],t=0..1,color=black,axes=normal);
Dk:=<diff(k(1),t),diff(k(2),t),diff(k(3),t)>: Dk:=simplify(Dk) assuming real:
DkNorm:=VectorNorm(Dk,2): DkNorm:=simplify(DkNorm) assuming real:
DkN:=Dk/DkNorm: DkN:=simplify(DkN) assuming real:
DDkN:=<diff(DkN(1),t),diff(DkN(2),t),diff(DkN(3),t)>: DDkN:=simplify(DDkN) assuming real:
DxB:=<diff(BReal(1),x),diff(BReal(2),x),diff(BReal(3),x)>: DxB:=simplify(DxB) assuming real:
DyB:=<diff(BReal(1),y),diff(BReal(2),y),diff(BReal(3),y)>: DyB:=simplify(DyB) assuming real:
DxBc:=subs(x=cReal(1),y=cReal(2),DxB): DxBc:=simplify(DxBc) assuming real:
DyBc:=subs(x=cReal(1),y=cReal(2),DyB): DyBc:=simplify(DyBc) assuming real:
N:=CrossProduct(DxBc,DyBc): N:=simplify(N) assuming real:
NN:=N/VectorNorm(N,2): NN:=simplify(NN) assuming real:
q:=Matrix([[1+NN(1)*I,NN(2)+NN(3)*I],[-NN(2)+NN(3)*I,1-NN(1)*I]]):
x:=Matrix([[DkN(1)*I,DkN(2)+DkN(3)*I],[-DkN(2)+DkN(3)*I,-DkN(1)*I]]):
vv:=MatrixMatrixMultiply(q,MatrixMatrixMultiply(x,MatrixInverse(q))):
v:=<Im(vv[1,1]),Re(vv[1,2]),Im(vv[1,2])>: v:=simplify(v) assuming real:
f:=v(1)*DDkN(1)+v(2)*DDkN(2)+v(3)*DDkN(3): f:=simplify(f) assuming real:
fS:=subs(t=s,f):
DkNormS:=subs(t=s,DkNorm):
sys:=diff(theta(s),s)=fS, diff(b1(s),s)=cos(theta(s))*DkNormS, diff(b2(s),s)=sin(theta(s))*DkNormS:
p:=dsolve({sys,theta(0)=0,b1(0)=0,b2(0)=0},{theta(s),b1(s),b2(s)},type=numeric):
plo:=odeplot(p,[b1(s),b2(s)],0..1,numpoints=200):
display(plo); |
Das Problem ist nun folgendes: Die Kurve c soll nun variieren. Nämlich will man verschiedene Winkel in der Definition von c oben. Im Code steht er momentan bei alpha=Pi/2. Wenn ich alpha=0 setze (also die Kurve ist dann c(t)=t), dann kommt bei der Berechnung von f weiter unten eine Art overflow. Wenn ich alpha=Pi/4 setze, dann dauert die Berechnung von DkNorm schon ewig. Sollte ich die Berechnung schon vor der Diffgleichung numerisch machen? Aber wie?
|
Profil
Quote
Link |
ThomasRichard
Aktiv
Dabei seit: 08.04.2010
Mitteilungen: 148
Aus: Aachen
 | Beitrag No.15, eingetragen 2012-09-14 11:35
|
Hallo calabi-yau,
diese Woche hat Oliver noch Urlaub, aber ich sehe, dass du auch auf MaplePrimes gefragt hast - gute Idee!
-----------------
Thomas Richard
Application Engineering / Technischer Support
Maplesoft Europe GmbH
|
Profil
Quote
Link |
calabi-yau
Senior
Dabei seit: 15.04.2006
Mitteilungen: 901
Aus:
| Beitrag No.16, vom Themenstarter, eingetragen 2012-09-15 09:59
|
Danke für die Info. Scheint eine schwierige Angelegenheit zu sein, auf MaplePrimes noch keine Antwort.
Ein kleines Zuckerl: Wenn ichs irgendwann hinbekomme, schreibe ich einen Artikel über das Modell auf dem MP. ;)
|
Profil
Quote
Link |
Ollie
Senior
Dabei seit: 03.05.2003
Mitteilungen: 5866
Aus: Aachen
| Beitrag No.17, eingetragen 2012-09-23 13:20
|
Hi,
Ich habe die Normroutine aus dem Algebrapaket nicht verwendet und einfach sqrt(M1^2+M2^2+M3^2) geschrieben. Das vermeidet overhead. Für manche Werte von alpha ergeben sich lange Rechenzeiten (z.B. Pi/5), wobei ich nicht weiß, ob das nicht einen mathematischen Grund (Nähe von Singularität) hat. Wenn keine Probleme auftreten, dauert ein Durchlauf bei mir etwa 25 Sekunden.
Gruß
Maple restart;ti:=time(): assume(x::real,y::real,t::real):
with(plots):with(LinearAlgebra):
M:=Re(1/(z^3-z^(-3)+sqrt(5))*<I*(z^2-z^(-2)),z^2+z^(-2),2*I/3*(z^3+z^(-3))>)+<0,0,2/5>:
B:=M/VectorNorm(M,2)^2:
BReal:=subs(z=x+I*y,B): BReal:=simplify(BReal) assuming real:
#plot3d([BReal(1),BReal(2),BReal(3)],x=-1..1,y=-1..1,grid=[80, 80]);
c:=exp(I*(3*Pi/4))*t;
cReal:=<Re(c),Im(c)>;
k:=subs(x=cReal(1),y=cReal(2),BReal): k:=simplify(k) assuming real:
#spacecurve([k(1),k(2),k(3)],t=0..1,color=black,axes=normal);
Dk:=<diff(k(1),t),diff(k(2),t),diff(k(3),t)>: Dk:=simplify(Dk,size) assuming real:
DkNorm:=(add(Dk(i)^2,i=1..3))^(1/2):DkNorm:=simplify(DkNorm) assuming real:
DkN:=Dk/DkNorm: DkN:=simplify(DkN) assuming real:
DDkN:=<diff(DkN(1),t),diff(DkN(2),t),diff(DkN(3),t)>:DDkN:=simplify(DDkN) assuming real:
printf("Zeit bis hierhin:%f \n",time()-ti);
DxB:=<diff(BReal(1),x),diff(BReal(2),x),diff(BReal(3),x)>: DxB:=simplify(DxB) assuming real:
DyB:=<diff(BReal(1),y),diff(BReal(2),y),diff(BReal(3),y)>: DyB:=simplify(DyB) assuming real:
DxBc:=subs(x=cReal(1),y=cReal(2),DxB): DxBc:=simplify(DxBc) assuming real:
DyBc:=subs(x=cReal(1),y=cReal(2),DyB): DyBc:=simplify(DyBc) assuming real:
N:=CrossProduct(DxBc,DyBc): N:=simplify(N) assuming real:
NN:=N/(add(N(i)^2,i=1..3))^(1/2): NN:=simplify(NN) assuming real:
q:=Matrix([[1+NN(1)*I,NN(2)+NN(3)*I],[-NN(2)+NN(3)*I,1-NN(1)*I]]):
x:=Matrix([[DkN(1)*I,DkN(2)+DkN(3)*I],[-DkN(2)+DkN(3)*I,-DkN(1)*I]]):
vv:=MatrixMatrixMultiply(q,MatrixMatrixMultiply(x,MatrixInverse(q))):simplify(vv):
vvc:=evalc(vv[1,2]):
v1:=Im(evalc(vv[1,1])):v2:=Re(vvc):v3:=Im(vvc):
v:=simplify~(v) assuming real:
f:=v1*DDkN(1)+v2*DDkN(2)+v3*DDkN(3): #f:=simplify(f) assuming real;
printf("Zeit bis hierhin:%f \n",time()-ti);
fS:=subs(t=s,f):
DkNormS:=subs(t=s,DkNorm):
sys:=diff(theta(s),s)=fS, diff(b1(s),s)=cos(theta(s))*DkNormS, diff(b2(s),s)=sin(theta(s))*DkNormS:
p:=dsolve({sys,theta(0)=0,b1(0)=0,b2(0)=0},{theta(s),b1(s),b2(s)},type=numeric):
plo:=odeplot(p,[b1(s),b2(s)],0..1,numpoints=200):
display(plo);
|
|
Profil
Quote
Link |
|
2013-05-25 22:06 U <
Wie unterstütze ich den Matheplaneten mit Geld?
[Neues Thema]