Snijpunten van Torussen

Wiskunde is niet alleen een vak op school. Kom je ergens in de praktijk (bijvoorbeeld tijdens je werk) een wiskundig probleem tegen dan kun je hier om hulp vragen.
Plaats reactie
Deuk
Nieuw lid
Nieuw lid
Berichten: 4
Lid geworden op: 07 sep 2015, 10:43

Snijpunten van Torussen

Bericht door Deuk » 07 sep 2015, 11:12

Ik ben op zoek naar een punt P in een 3d ruimte. Het enige hulpstuk wat ik heb is een driehoek op de 'grond' met de hoeken ABC. Van deze driehoek zijn alle coordinaten bekend en van punt P zijn alleen de hoeken met driehoek ABC bekend (dus hoek APB, APC, en BPC).
Afbeelding

In 2D en met de regel van de ingeschreven hoek kan punt P op ieder punt van een boog door de punten van een ribbe van de driehoek liggen, zie bijvoorbeeld DPC in dit plaatje van Wikipedia:
Afbeelding

In 3D betekend dat dus dat je torussen (wiskundige donuts, in dit geval zonder gat) kan draaien om de 3 ribben, en dat punt P dus op ieder snijpunt van die drie torussen kan liggen. Maar nu is de vraag: hoe vind ik de snijpunten van torussen?

Met CAD software heb ik kunnen "bewijzen" dat het punt P inderdaad op een van de vier snijpunten ligt. Maar wiskundig gezien kom ik er maar niet uit. Mocht ik eenmaal de snijpunten gevonden hebben dan is het wegstrepen van de 'foute' snijpunten geen probleem aangezien het punt P zich nooit in de negatieve ruimte bevindt.

Afbeelding Afbeelding

Ik heb de vraag op math.stackexchange gepost, maar daar is er nog niet op gereageerd. Voor meer uitleg en een poging om met de formules van torussen aan de slag te gaan zie:
http://math.stackexchange.com/questions ... ri-toruses

Mocht je me kunnen helpen hoor ik dat graag!
Bij voorbaat dank.

Groet,
Paul

arie
Moderator
Moderator
Berichten: 3911
Lid geworden op: 09 mei 2008, 09:19

Re: Snijpunten van Torussen

Bericht door arie » 07 sep 2015, 20:22

Als je een numerieke oplossing wil, dan kan dat via de cosinusregel:
noem:
AB = a
BC = b
CA = c
AP = d
BP = e
CP = f
hoek APB = alpha
hoek BPC = beta
hoek CPA = gamma

Uit de gegeven coordinaten van A, B en C bereken je a, b en c.
Ook kennen we alpha, beta en gamma.
De waarden van d, e en f willen we weten, immers:
P ligt op de bollen:
[1] met middelpunt A en straal d
[2] met middelpunt B en straal e
[3] met middelpunt C en straal f
en 3 bollen kunnen we relatief eenvoudig snijden.

Via de cosinusregel gelden de volgende 3 gelijkheden:








Uit de middelste van deze 3 volgt:



dus via de abc-formule:



ofwel




Evenzo uit de laatste van de 3 cosinusregel vergelijkingen:




Herschrijven we de eerste van de 3 cosinusregel vergelijkingen:



en substitueren we hierin de bovenstaande resultaten voor d en e (2 mogelijkheden per variabele), dan houden we links een uitdrukking over waarin alleen f onbekend is. Het nulpunt hiervan is numeriek te bepalen, waarna we de waarde van f kennen.
En we hadden d en e beide al als functie van f, dus dan kennen we die 2 ook.
Resteert nog het snijden van de 3 bollen.

Kom je hiermee verder?

Deuk
Nieuw lid
Nieuw lid
Berichten: 4
Lid geworden op: 07 sep 2015, 10:43

Re: Snijpunten van Torussen

Bericht door Deuk » 08 sep 2015, 14:25

Pfoe, ik ben bang van niet! :(
Het substitueren van d en e moet nog wel lukken, maar hoe ik dan f daarmee moet vinden zou ik zo niet weten.

Zelf heb ik een andere aanpak geprobeerd door niet de lengtes d, e en f als radius van een bol te gebruiken maar door bogen om de ribben a, b en c te wentelen. Dan krijg je de torussen in 3D ruimte. Ik ben inmiddels zo ver dat ik de kleine en grote radius van de torus kan berekenen:
Afbeelding

Met de definitie van een torus zou ik dan moeten kunnen berekenen wat alle mogelijk posities van x, y en z zijn voor het punt p in de drie torussen. Met de translatie in 3D space kom ik volgens mij uit op deze formules, waarbij opposite de lengte is van ribbe a, b en c, uitgaande van een gelijkbenige driehoek:

Afbeelding













Nu weet ik alleen niet 100% zeker of deze formules kloppen, maar als ik op basis van deze formules torussen genereer in matlab lijkt dat wel zo te zijn. Ik heb hier als voorbeeld punt p op 40, 70, 100 gezet, en met meten in 3D programma Rhino geeft dat aan dat de hoeken 22.83, 23.81 en 21.85 zijn. De complete matlab code is hier terug te vinden: http://math.stackexchange.com/questions ... ri-toruses
Afbeelding

Mijn idee was om de drie snijpunten te vinden door de vergelijkingen x1 = x2 = x3, y1 = y2 = y3, en z1 = z2 = z3 op te lossen. Maar dat ging mijn wiskundige petje te boven. Ik had een stukje code geschreven die automatisch door alle waarden van l en k heen liep om zo te kijken waar ze snijden maar dat duurde enooooorm lang (uren) en ik kon volgens mij wel de waarden van x en y vinden (40, 70), maar de waarde van z paste daar dan weer niet bij (80 of iets dergelijks?).

Klinkt dit als een goede aanpak en zou dit opgelost kunnen worden?
Of is het verstandiger om f te vinden en dan de parametrische vergelijking van een bol te gebruiken? (maar dan heb ik het zelfde probleem als wat ik nu heb met de torussen toch?)

arie
Moderator
Moderator
Berichten: 3911
Lid geworden op: 09 mei 2008, 09:19

Re: Snijpunten van Torussen

Bericht door arie » 10 sep 2015, 08:57

Je bent uitgegaan van 2 parameters, l en k.
Echter: elke torus heeft zijn eigen parameters, dus:
- torus 1: l1 en k1
- torus 2: l2 en k2
- torus 3: l3 en k3
Immers: het is niet zo dat als de parameters van 2 torussen gelijk zijn, dat dan ook het punt (x,y,z) op beide torussen gelijk is. Anders zou dat voor elk tweetal (l,k) moeten gelden, en zijn ook de 2 torussen hetzelfde.

Je code met aangepaste parameters wordt dus:

Code: Selecteer alles

x1 = (sind(90)*minorR1*sind(l1)) + opposite1 / 2;
y1 = (majorR1+minorR1*cosd(l1)).*sind(k1);
z1 = -sind(90)*(majorR1+minorR1*cosd(l1)).*cosd(k1);

x2 = (sind(90) * (minorR2*sind(l2))) * cosd(120) - ((majorR2+minorR2*cosd(l2)).*sind(k2)) * sind(120) + sind(30)*opposite/2 + opposite/2;
y2 = (sind(90) * (minorR2*sind(l2))) * sind(120) + ((majorR2+minorR2*cosd(l2)).*sind(k2)) * cosd(120) + (sind(60)*opposite / 2);
z2 = -sind(90) * ((majorR2+minorR2*cosd(l2)).*cosd(k2));

x3 = (sind(90) * minorR3*sind(l3))*cosd(240) - ((majorR3+minorR3*cosd(l3)).*sind(k3)) * sind(240) + (sind(30) * opposite3 / 2);
y3 = (sind(90) * minorR3*sind(l3))*sind(240) + ((majorR3+minorR3*cosd(l3)).*sind(k3)) * cosd(240) + (sind(60) * opposite3 / 2);
z3 = -sind(90) * (majorR3+minorR3*cosd(l3)).*cosd(k3);

Als verder x1=x2=x3=x, y1=y2=y3=y en z1=z2=z3=z hou je 9 vergelijkingen over met 9 onbekenden (x, y, z, k1, l1, k2, l2, k3, l3).
Dit stelsel is niet lineair: het bevat producten van variabelen en bovendien ook sinus- en cosinusvormen.
Ik verwacht daarom dat oplossen van dit stelsel niet eenvoudig zal zijn.

Om 6 parameters (l1, l2, l3, k1, k2, k3) volledig af te zoeken kost nog veel meer tijd dan je 2 parameters.
Over het algemeen is uitputtend zoeken zeer bewerkelijk, en zal je snellere zoekmethoden willen gebruiken, bijvoorbeeld:
[1] begin met een punt P0 = (xp, yp, zp), bepaal de totale fout tussen elk van de 3 gewenste hoeken (APB, BPC, CPA) en de huidige hoeken (AP0B, BP0C, CP0A)
[2] kijk vervolgens naar een aantal buurpunten in de directe omgeving, bijvoorbeeld

en bepaal van elk van die punten hun fout t.o.v. de gewenste hoeken.
[3] neem als volgende punt Pn het buurpunt met de kleinste fout.
[4] zolang de fout afneemt: ga door via punt [2]



Als alternatief voor je stelsel heb ik hierboven d en e als functie van f geschreven:





Substitutie in



levert een vergelijking met f als enige onbekende.
NOOT: dit is niet de enige en/of beste methode, maar volgens mij wat eenvoudiger om op te lossen.

In Pari/GP is de code hiervoor:

Code: Selecteer alles

\\ Euclidean norm of vector v:
eunorm(v)={
sqrt(v*v~)
}

\\ cross product of the 3D vectors u and v:
cross(u,v) = {
s = vector(3);

s[1] = u[2]*v[3] - u[3]*v[2];
s[2] = u[3]*v[1] - u[1]*v[3];
s[3] = u[1]*v[2] - u[2]*v[1];
s
}

\\ return d as function of f:
fnd(f)={
f*cos(gam) + sqrt(c^2 - (f*sin(gam))^2)
}

\\ return e as function of f:
fne(f)={
f*cos(bet) + sqrt(b^2 - (f*sin(bet))^2)
}


\\ MAIN PROGRAM:
{
\\ given values:
A = [0,0,0];
B = [50,0,0];
C = [25, 25*sqrt(3), 0];
alpd = 22.82699023521976221689472299;
betd = 23.80770096561971797946120879;
gamd = 21.85078089394426844695048013;

\\ calculate the base triangle side lengths:
a = eunorm(A-B);
b = eunorm(B-C);
c = eunorm(C-A);

\\ convert the given angles to radians:
alp = alpd*Pi/180;
bet = betd*Pi/180;
gam = gamd*Pi/180;

\\ calculate upper bound of f:
fmaxb = b/sin(bet);
fmaxc = c/sin(gam);
fmax = min(fmaxb, fmaxc)*0.999;

\\ find zero (numerical solution):
f = solve(f = 0.001, fmax, fnd(f)^2 + fne(f)^2 - a^2 - 2*fnd(f)*fne(f)*cos(alp));

\\ now we have f, calculate d and e:
d = fnd(f);
e = fne(f);

\\ show radii:
print("radii:");
print("d = ", d);
print("e = ", e);
print("f = ", f);


\\ intersect the spheres:
vBA = B-A;
vCA = C-A;
cp = cross(vBA, vCA);
c2 = cp*cp~;
fCA = vBA*vBA~ + d^2 - e^2;
fBA = vCA*vCA~ + d^2 - f^2;
u = cross( (fCA*vCA - fBA*vBA)/2, cp) / c2;
v = sqrt(d^2-u*u~)*cp/sqrt(c2);
sol1 = A+u+v;
sol2 = A+u-v;

print("\nSolutions:");
print("P1 = ", sol1);
print("P2 = ", sol2);
}
Eerst zie je 2 hulpfuncties: eunorm() en cross(), volgens mij zitten die in Matlab als norm() en cross().
Dan
fnd(f) = d als functie van f
fne(f) = e als functie van f
en vervolgens het hoofdprogramma:
- eerst de gegeven (=bekende) waarden
- dan wat waarden van zijden en hoeken die we daaruit direct kunnen afleiden
- dan de bovengrens van f (in bovenstaande formules willen we geen imaginaire getallen: onder het wortelteken moet niet-negatief zijn)
- dan de solve() functie: deze functie vindt nulpunten van een gegeven (continue) functie, dus een oplossing van
fun(x) = 0,
waarbij variabele x in een interval [a, b] ligt.
De code om een nulpunt x0 te vinden wordt dan:
x0 = solve( x = a, b, fun(x) )
ENIGE VOORWAARDE: fun(a) en fun(b) moeten aan verschikllende kanten van de x-as liggen.
Volgens mij heeft Matlab ook zoiets, zoek daarin eens bij solve, fsolve, fzero of fnzero.
Zo niet, dan is het niet al te ingewikkeld om zelf zo'n functie te maken.
- uit f volgen ook d en e
- het snijden van 3 bollen is min of meer standaard-code

In Pari/GP is
u * v het inproduct van vectoren u en v, en
v~ de getransponeerde van v

De output van dit programma is:

radii:
d = 128.4523257866512902011689848
e = 122.4744871391589049098642038
f = 104.5840435893971804344229717

Solutions:
P1 = [40.00000000000000000000000000, 69.99999999999999999999999989, 100.00000000
00000000000000001]
P2 = [40.00000000000000000000000000, 69.99999999999999999999999989, -100.0000000
000000000000000001]


NOOT: fnd() en fne() kijken elk naar 1 van hun 2 mogelijke oplossing, de andere 3 combinaties moet je ook nog evalueren.

Deuk
Nieuw lid
Nieuw lid
Berichten: 4
Lid geworden op: 07 sep 2015, 10:43

Re: Snijpunten van Torussen

Bericht door Deuk » 15 sep 2015, 10:25

Je laat het zo makkelijk lijken, ∞ thanks voor je hulp!!

Ik heb geprobeerd om het naar Matlab te vertalen, maar dat ging niet echt. Toen gekeken naar Pari/GP en dat geïnstalleerd, dat kreeg ik wel aan de praat, werkt perfect!

Over het evalueren van die andere oplossingen, is het daarvoor een idee om twee "kopieën" van fnd en fne te maken waarbij de uitkomst van de wortel negatief is, en dan 3 solve functies toe te voegen met de drie mogelijke opties? En pas als dat dan een valide oplossing oplevert voor f de rest van de formules uitvoeren. Of is dat te kort door de bocht?

Code: Selecteer alles

\\ Euclidean norm of vector v:
eunorm(v)={
sqrt(v*v~)
}

\\ cross product of the 3D vectors u and v:
cross(u,v) = {
s = vector(3);

s[1] = u[2]*v[3] - u[3]*v[2];
s[2] = u[3]*v[1] - u[1]*v[3];
s[3] = u[1]*v[2] - u[2]*v[1];
s
}

\\ return d as function of f:
fnd(f)={
f*cos(gam) + sqrt(c^2 - (f*sin(gam))^2)
}

fndneg(f)={
f*cos(gam) - sqrt(c^2 - (f*sin(gam))^2)
}

\\ return e as function of f:
fne(f)={
f*cos(bet) + sqrt(b^2 - (f*sin(bet))^2)
}

fneneg(f)={
f*cos(bet) - sqrt(b^2 - (f*sin(bet))^2)
}

sphereintersect(f)={
	\\ now we have f, calculate d and e:
	d = fnd(f);
	e = fne(f);

	\\ show radii:
	print("radii:");
	print("d = ", d);
	print("e = ", e);
	print("f = ", f);


	\\ intersect the spheres:
	vBA = B-A;
	vCA = C-A;
	cp = cross(vBA, vCA);
	c2 = cp*cp~;
	fCA = vBA*vBA~ + d^2 - e^2;
	fBA = vCA*vCA~ + d^2 - f^2;
	u = cross( (fCA*vCA - fBA*vBA)/2, cp) / c2;
	v = sqrt(d^2-u*u~)*cp/sqrt(c2);
	sol1 = A+u+v;
	sol2 = A+u-v;

	print("\nSolutions:");
	print("P1 = ", sol1);
	print("P2 = ", sol2);
	print("\n\n\n");
}


\\ MAIN PROGRAM:
{
	\\ given values:
	A = [0,0,0];
	B = [50,0,0];
	C = [25, 25*sqrt(3), 0];

	alpd = 22.82699023521976221689472299;
	betd = 23.80770096561971797946120879;
	gamd = 21.85078089394426844695048013;

	\\ calculate the base triangle side lengths:
	a = eunorm(A-B);
	b = eunorm(B-C);
	c = eunorm(C-A);

	\\ convert the given angles to radians:
	alp = alpd*Pi/180;
	bet = betd*Pi/180;
	gam = gamd*Pi/180;

	\\ calculate upper bound of f:
	fmaxb = b/sin(bet);
	fmaxc = c/sin(gam);
	fmax = min(fmaxb, fmaxc)*0.999;

	\\ find zero (numerical solution):
	f = solve(f = 0.001, fmax, fnd(f)^2 + fne(f)^2 - a^2 - 2*fnd(f)*fne(f)*cos(alp));
	sphereintersect(f);

	f = solve(f = 0.001, fmax, fndneg(f)^2 + fne(f)^2 - a^2 - 2*fndneg(f)*fne(f)*cos(alp));
	sphereintersect(f);

	f = solve(f = 0.001, fmax, fnd(f)^2 + fneneg(f)^2 - a^2 - 2*fnd(f)*fneneg(f)*cos(alp));
	sphereintersect(f);

	f = solve(f = 0.001, fmax, fndneg(f)^2 + fneneg(f)^2 - a^2 - 2*fndneg(f)*fneneg(f)*cos(alp));
	sphereintersect(f);


}
Wel moet ik hier nog afvangen of de solve geen domain error geeft. Dat ga ik nog even uitzoeken.

arie
Moderator
Moderator
Berichten: 3911
Lid geworden op: 09 mei 2008, 09:19

Re: Snijpunten van Torussen

Bericht door arie » 15 sep 2015, 20:29

Mooi dat het werkt.

Die foutmelding wordt waarschijnlijk veroorzaakt doordat in sommige gevallen je domeingrenzen een functiewaarde met hetzelfde teken geven.
Hierop kan je testen via de sign() functie van Pari/GP.

Wel vreemd dat het in Matlab niet lukt, we gebruiken geen heel bijzondere rekenkundige functies.

Deuk
Nieuw lid
Nieuw lid
Berichten: 4
Lid geworden op: 07 sep 2015, 10:43

Re: Snijpunten van Torussen

Bericht door Deuk » 16 sep 2015, 09:24

Dat het in Matlab niet werkt ligt meer aan mijn ervaring met Matlab dan de functionaliteit van Matlab!

Wederom bedankt! :D

Plaats reactie