Konstrukcija elipse

[Medo5]

Super, hvala.

Čisto sem pozabil na piko .

Za rezultat matrike sem pričakoval podobne številke, kot jih dobim z mojo prejšno. Pri meni je več samo to, da dobim še zadnji člen in potem bom vse prejšne člene delil s to coifro krat (-1).
Nadalje sem številke vstavil v tvojo matriko B in vrne pravilen razultat .

Vse, kar moram še storiti je poiskati lastne vrednosti matrike B. Glede nato, da je 2x2 upam, da obstaja kar končna enačba in bom hitro spisal.

Ko končam pošljem rezultat.

Lp; Jan

[Aniviller]

Seveda obstaja koncna enacba - nekaj sem ti jih ze izpeljal, problem je le, da obstaja par posebnih primerov pri katerih bos dobil 0/0. Te bos rocno preverjati. Ti primeri so najverjetneje tisti, pri katerih je elipsa ze poravnana. Enako velja za iskanje sredisca, tam so tudi posebni primeri.
Pa pazi na to, da prosti clen (F) se vpliva na velikost elipse.

Lastne vektorje tudi dobis z direktno formulo (spet je posebni primer ko je elipsa ze v smeri osi). Potem samo tistega od vecje lastne vrednosti normiras in mu poisces kot.

[Medo5]

Vse sem že skodiral, edino kot mi še vedno dela preglavice. Pri splošni enačbi sta še vedno problematiča člena a in c, ki sta enaka.
Je možno dobiti kot že iz matrike B?

Lp; Jan

[Aniviller]

kot ne bi smel biti problem. c je misljen clen pri xy, a in b pa pripadata kvadratom. Tudi v matriki B je izvendiagonalni element seveda tisti od clena xy. Mislim da sem res malo spreminjal oznake. Pa tudi tista formula ni cisto prav, sem hitro pisal. Ce imas
\(\phi=\arctan\frac{2(a-\lambda)}{c}\)
je vse ok. (razen ce je c 0, takrat je itak tudi kot 0). \(\lambda\) je ena lastna vrednost (in ne polos, kot sem napisal prej). Ekvivalentna je formula
\(\phi=\arctan\frac{c}{2(b-\lambda)}\)
Odvisno pac katero os vzames za glavno in kako imas ostevilcene lastne vrednosti. Najvec kar se ti lahko zgodi je, da moras kotu na koncu odsteti ali pristeti pravi kot

[Medo5]

Malo z zamudo, ampak vseeno pošiljam rezultate.
Prva slika: izračun elipse na podlagi petih točk, ki se vzamejo iz polinoma. Problem je, ker se polinom ne uspe lepo napeljat po točkah dobljenih iz te rotirane elipse. Polinom je stopnje 10, več negre, ker ga uniči numerična napaka. Če je elipsa bolj podobna krogu, se polinom lepo napelje in tudi rezultati so vredu. Tako sliko sem dal samo za primerjavo glede na drugi algoritem.

Druga slika: elipsa se izračuna iz vseh 360 točk, po metodi najmanjših kvadratov, ki jo je nakazal Aniviller. Razlika je očitna.

Šenkrat najlepša hvaza za pomoč. Zanima pa me še, od kej ti vso to matematično znanje in predvsem razumevanje vsega skupaj?

Lp; Jan

Iz petih točk na polinomu

Iz vseh točk

[Aniviller]

Bravo, tole ti je pa lepo uspelo. Sem se bal da bodo kaksni problemi s predelavo podatkov.

Studiram fiziko, matematicna smer, kar je natanko to, kar rabis za tvoj problem Studij fizike daje najvec matematicnega znanja v primerjavi z ostalimi fakultetami (pustimo same matematike ob strani), in tudi vsega uporabimo na fizikalnih problemih, kjer je potrebno razumevanje. Tako abstraktne matematicne postopke lahko povezes s prakticnim problemom (najbolj vektorski in matricni racun, lastni sistemi, razvoj v vrste (se pravi celotna linearna algebra), integrali, diferencialne enacbe, minimizacija, v matematicni smeri pa poglobljene numericne metode za vse nasteto).

eno vprasanje, del kaksnega programa pa je tole? izgleda zanimivo.

[NIKKI]

Bravo!
Res si !

[Medo5]

Najprej hvala obema za pohvale .

Sam študiram na FE, ampak žal od tega faksa nebom odnesel nič. Preveč teorije, nobenemu se nič neda, nič uporabnega v praksi. No, mogoče sedaj ko sem že blizu konca, ampak sedaj je že prepozno.
Za nadalnji študi pa resno razmišljam o FMF ali pa FRI.

Končni izdelek je precizno merjenje premerov nekih izvrtin. Priložene slike so iz Corela, zato so "popolne". V resnici je na sliki veliko motenj, res pa je, da so odstopnja od koga majhna, zato jih lahko izločim.

Imam pa šeeno vprašanje:
Pri drugi metodi je zahtevan PseudoInverse. Kakor sem si uspel pogledat obstaja nekaj načinov izračuna (SVD,QR,..).
Nazadnje sem enostavno naredil po definiciji pseudoinverza:B^+=(B^(T)B)^(-1)B^(T).
Pomanjklivost je ta, da če je matrika singularna, nima inverza. Zanima me, kakšna ja verjetnost, da je matrika s 100-timi točkami tipa double, singularna? Upam, da ne se motim, ampak mislim, da se v praksi to ne more zgoditi .

Kakšne so bistvene prednosti SVD,QR?

Lp; Jan

[Aniviller]

Tvoja metoda je za ta namen najboljsa - na lepse se napise to kot resevanje matricne enacbe:
\(Bx=b\to B^T Bx=B^T b\)
kar je seveda isto.
SVD razcep je bistveno splosnejsi in da cel kup dodatnih podatkov, ki jih sploh ne potrebujes. Pomen SVD-ja je, da ti matriko razcepi na vsoto matrik, ki se lahko zapisejo kot produkt dveh vektorjev. S tem razcepis matriko na pomembne in manj pomembne dele.
QR razcep ti ne prinese nic novega: itak za samo racunanje inverza matrike potrebujes QR ali LU razcep (to tudi interno uporabljajo vse metode obracanja matrik). Nihce namrec ne obraca matrik po definiciji.

Mislim da je opisana metoda kar najboljsa. Tudi ko bos preklopil na grse podatke, mislim da ne bo tezav, ker metoda najmanjsih kvadratov ni tako obcutljiva na slabe podatke. Mislim, da lahko tudi zmanjsas stevilo tock po kotu (ceprav te nic ne kosta ce jih je vec, ker matrika, ki jo obracas ni nic vecja).

Ce je pa matrika \(B^T B\) singularna, pa pomeni da problem itak nima prave resitve. To se ti ne bo zgodilo (ce pa se, se pa ne da pomagat).

Glede studija pa lahko recem, da na FRI ne bos izvedel nic novega. Ce isces programiranje, potem je to matematicna smer na fiziki. Manjkajo le bolj strukturno naravnani algoritmi (drevesa in drugo ''rastlinje'') in stroga racunalniska teorija, ki pa je itak ena najbolj dolgocasnih tem (ce te pa zanima, pa tudi ni tako komplicirano, da si ne bi mogel sam pogledat). Moram pa opozoriti, da je zacetni preboj precej tezek, pa tudi kasneje ni hec.

[Medo5]

Ker je tip problema podoben obstoječemu, bom nadaljeval kar v tej temi.

Zadnje dni se ukvarjam s fitanjem cilindra na množico točk v 3D prostoru.
Izhodiščna enačba je splošna enačba "quadric surface" .

Koda: Izberi vse

Ax2 + By2 + Cz2 + Dxy + Exz + Fyz  + Gx + Hy + Iz +  J = 0

Enačba ima 10 členov, od tega lahko damo konstanto na npr 1 in s tem imamo končno število parametrov 9.
Osnovna ečaba eliptičnega cilindra je

Koda: Izberi vse

x2/a2 + y2/b2 =1 

Ima samo dve neznanki.
Nikjer ne najdem enačbe, ki bi vsebovala vseh devet neznank.
Možne neznanke ki se jih spomnim so: X0, Y0, Z0, r1, r2, alfa, beta, gama. Ena mi še vedno manjka.
S pomočjo pseudo inverza, ki ga je opisal Aniviller, lahko dobim koeficiente A,B,C,...
Vendar zaenkrat nevem, kako bi jih povezal. Predvsem me moti, ker bi imel več enačb, kot pa neznank.
V Mathematici sem poizkusil zgenerirati osnovne enačbe, kjer sem upošteval X0, Y0, gama, beta. Torej zamaknjen center in rotacijo v prostoru po osi X in Y. Ker je cel kup trigonometričnih funkcij, se mi zdi sistem nerešljiv. Sploh pa, ker imam preveč enačb.

Koda: Izberi vse

Rxg = {
   {1, 0, 0},
   {0, Cos[Gamma], -Sin[Gamma]},
   {0, Sin[Gamma], Cos[Gamma]}

Ryb = ...
Rza = ...

Baza = {{x},{y},{z}}

Simplify[Rxg.Ryb.Baza] // MatrixForm
xx = x Cos[\[Beta]] + z Sin[\[Beta]]
yy = y Cos[\[Gamma]] + (-z Cos[\[Beta]] + x Sin[\[Beta]]) Sin[\[Gamma]]

Expand[(xx - X0)^2 + (yy - Y0)^2  - r^2 ]
-r^2 + X0^2 + Y0^2 - 2 x X0 Cos[\[Beta]] + x^2 Cos[\[Beta]]^2 - 
 2 y Y0 Cos[\[Gamma]] + y^2 Cos[\[Gamma]]^2 - 2 X0 z Sin[\[Beta]] + 
 2 x z Cos[\[Beta]] Sin[\[Beta]] + z^2 Sin[\[Beta]]^2 + 
 2 Y0 z Cos[\[Beta]] Sin[\[Gamma]] - 
 2 y z Cos[\[Beta]] Cos[\[Gamma]] Sin[\[Gamma]] - 
 2 x Y0 Sin[\[Beta]] Sin[\[Gamma]] + 
 2 x y Cos[\[Gamma]] Sin[\[Beta]] Sin[\[Gamma]] + 
 z^2 Cos[\[Beta]]^2 Sin[\[Gamma]]^2 - 
 2 x z Cos[\[Beta]] Sin[\[Beta]] Sin[\[Gamma]]^2 + 
 x^2 Sin[\[Beta]]^2 Sin[\[Gamma]]^2
  

V primeru elipse je bilo 5 enačb, 5 neznak, zato se je vse skupaj izšlo.

Vsak namig bo zelo dobrodošel, še posebej, ker sem mnenja, da sem verjetno nekaj generalnega pozabil. Ali pa so moje želje prevelike in preoptimistične.

Lp; Jan

[Aniviller]

Kot v prejsnjem primeru, tukaj tudi nafitaj na prvo enacbo... potem pa pride do iskanja lastnih vrednosti in lastnih vektorjev 4x4 matrike (4 zato, ker je se premaknjeno). Dobis dve pozitivni lastni vrednosti za lastne smeri elipse in eno nicelno za v smeri cilindra (ta ne bo tocno nic zaradi napak - ta bo lahko sluzila kot test, kako pravilno cilindricno je v resnici). Ena lastna vrednost bo pa 1 (ta bo pripadala premikom).

Trenutno nimam casa, ti bom v polnosti odgovoril enkrat cez dan.

[Aniviller]

Uvedba kotov je slaba ideja... pri kotih takoj dobis notri kotne funkcije, pa prednost das dolocenim smerem. Najboljsa resitev bo kar fitanje na splosno funkcijo, kot si napisal (ta postopek za dobit A,B,C,... je itak isti kot prej). Potem moras pa dolocit sredisce te ploskve, da lahko premaknes v izhodisce, ter dolocis lastne smeri in polosi (tisto kar sem prej govoril o 4x4 matriki ni dobro). Problem: ce je res cilinder, sredisce ni dobro doloceno (kjerkoli na osi je lahko). Zato bo treba uporabljati psevdoinverze, da dobis pac neko resitev (vse so enako dobre).

Torej, definirajmo matriko M:
\(M=\begin{bmatrix}A&D/2&E/2\\D/2 &B&F/2\\E/2&F/2&C\end{bmatrix}\)
Sredisce dobis iz sistema enacb
\(M\begin{bmatrix}x_0\\y_0\\z_0\end{bmatrix}=-\begin{bmatrix}G/2\\H/2\\I/2\end{bmatrix}\)
Predlagam PseudoInverse[] od Mathematice. S tem imas sredisce.

Ce zamenjamo \(x=x'+x_0\) enacba postane
\(A{x'}^2+B{y'}^2+C{z'}^2+Dx'y'+Ex'z'+Fy'z'+\)\((J+Ax_0^2+By_0^2+Cz_0^2+Dx_0y_0+Ex_0z_0+Fy_0z_0)=0\)

Kot vidimo, tip ploskve, ki ga dolocajo A do F, ni odvisen od tistih linearnih clenov, vplivajo pa na debelino (vse tisto v oklepaju je zdaj prosti clen). Preimenujmo ta oklepaj:
\(Ax^2+By^2+Cz^2+Dxy+Exz+Fyz=Q\)

Lastne vrednosti matrike M so povezane s polosmi kot zadnjic (ce je cilinder, bo ena morala biti nic... oziroma skoraj nic - neko odstopanje bo, ampak to je v redu ker imas s tem dober test kako cilindricna je stvar v resnici). Drugi dve lastni vrednosti bosta pa predstavljali polosi elipse v preseku. Za dejanske velikosti polosi bo prisel v igro tisti prosti clen Q (\(a=\sqrt{Q/\lambda}\) za vsako polos a s pripadajoco lastno vrednostjo \(\lambda\)).

Pripadajoci lastni vektorji so pa smeri - tisti od nicelne lastne vrednosti bo smer osi cilindra.

[Medo5]

Hvala za odgovor. Nekaj časa bom potreboval, da si popolnoma razjasnim stvari, do takrat pa prilagam rezultat dobljen v Mathematici. Je točno to, kar iščem, vednar je samo prikaz. Opazil sem tudi, da je zelo majhna verjetnost, da res pride eliptičen cilinder. Nisem se še poglabljal, vendar ali je možno enačbo nastaviti tako, da bo rezultat vedno eliptičen cilinder, pa naj bodo vhodne točke karkoli. Seveda ujem nebo popolen.

Koda: Izberi vse

PM = PseudoInverse[A];
V = Table[1, {12}];
res = N[PM.V]
eq = x^2*res[[1]] + y^2*res[[2]] + z^2*res[[3]] + x*y*res[[4]] + 
   x*z*res[[5]] + y*z*res[[6]] + x*res[[7]] + y*res[[8]] + 
   z*res[[9]] - 1;

eqn = eq == 0;
s = ContourPlot3D[eq == 0, {x, 0, 10}, {y, 0, 10}, {z, -20, 20}];
p = ListPointPlot3D[pts, PlotStyle -> {PointSize[0.05]}];
Show[s, p]

Ali bobstaja enačba cilindra, ki bi vsebovala vse te parametre, ki jih bom dobil iz koeficientov A,B,C,..?

Lp; Jan

[Aniviller]

Seveda je neskoncno majhna verjetnost da res dobis cilinder (ker moras od vseh realnih stevil okrog nicle tocno niclo zadet in zato dobis zelo dolg elipsoid ali zelo dolg hiperboloid - gre malo narazen ali skupaj). Vendar to ni problem: polosi dobljene elipse sta vseeno pravilni in tisto lastno vrednost lahko vedno nastavis nazaj na nic.

Res bo bolj pregledno ce s fitom dobis tiste A...J, dobis sredisce, potem pa res naredis matriko M in poisces lastne vrednosti in lastne vektorje. Ce matrika M predstavlja tvojo enacbo (sicer v teziscnem sistemu), in ce jo zapises kot \(M=PDP^{-1}\) (P je matrika lastnih vektorjev (to je tisto kar ti da ven Eigensystem), D pa diagonalna matrika z lastnimi vrednostmi po diagonali). Ce zdaj temu D-ju "skoraj niclo" spremenis v niclo, dobis matriko z novimi koeficienti A...F, ki predstavljajo pravi cilinder, ki se najbolje prilega tvojim podatkom. Po domace povedano, najdes lastni koordinatni sistem, zravnas v cilinder in obrnes nazaj.
S tem postopkom bos vedno dobil koeficiente A do F, ki predstavljajo cilinder.

[Medo5]

Sem uspel sestavit skupaj. V enačbah sem našel nekaj malega napak, bom na koncu zbral skupaj, če se bo še kdo v prihodnosti ukvarjal s tem.
Trenutno me daje še problem s končnim fitanjem cilindra. Kot si napisal v realnosti ni mogoče dobiti ravno cilindra. Zato sem v postopek dodal še tvoj zadnji namig, kako dobiti cilinder. Rezultat je res cilinder, vendar ni nujno da je "pravi" cilinder. Trenutno dobim Hyperbolic cylinder, ki pa ima drugačne vrednosti za (center, radij,..) kot pa moj željeni cilinder. Ali je možnost, da dobim ravno eliptični cilinder?

Na splošno pa bi ponovno rad pohvalil tvoje super matematično znanje in se zahvalil za pomoč .

Lp; Jan