Konstrukcija elipse

O matematiki, številih, množicah in računih...
Medo5
Posts: 26
Joined: 8.2.2007 9:08
Location: Ljubljana
Contact:

Re: Konstrukcija elipse

Post by Medo5 » 7.1.2010 12:38

Najprej bom izpostavil moje spremembe in upam, da sem postopal pravilno.
V enačbo za Q sem namesto x,y,z vstavil dobljeno središče x0, y0, z0 in prištel konstanto 1. Torej Q= A*x0^2 + ....+F*y0*z0 + 1

Za računanje polosi sem vzel absolutno vrednost pod korenom. V nasprotnem primeru lahko pride do korenjenja negativnega števila.

Code: Select all

A->vhodna matrika
PM = PseudoInverse[A];
V = Table[1, {48}];
res = PM.V
M = ....
MI = PseudoInverse[M];
B={res[[7]]/2,res[[8]]/2,res[[9]]/2}
CENTER = -MI.B;

R1 = Sqrt[Abs[Q/lambda[[1]]]]
R2 = Sqrt[Abs[Q/lambda[[2]]]]
Sedaj pa še zadnja iteracija:

Code: Select all

lambda = Eigenvalues[M]
 vector = Eigenvectors[M]
 P={{vector [[1, 1]], vector [[1, 2]], vector [[1, 3]]}, {vector [[2, 1]], vector [[2, 2]], vector [[2,3]]},
{vector [[3, 1]], vector [[3, 2]], vector [[3, 3]]}}
PInv = Inverse[P];
Diag ={{lambda[[1]], 0, 0}, {0, lambda[[2]], 0}, {0,0,0}}
M2 = PInv.Diag.P;
MM = M2.CENTER*-2;
Sprememba je v izračunu M2. Pomešan je vrstni red členov. Delal sem po principu, da če vstavim notri vse tri lambde, potem moram dobiti enak rezultat, kot pred zadnjo iteracijo.

Rezultat tega je hiperbolični cilinder(za moj primer točk).
Na spletu sem našel seznam, kako prepoznati za katero obliko gre; link: http://mathworld.wolfram.com/QuadraticSurface.html .
Iz tabele je razvidno, da ima eliptični cilinder obe lastni vrednosti enakega predznaka. Če v zadnji iteraciji upoštevam še to, potem je rezultat že zelo blizu željenega. Ali pa sem imel samo precej sreče. Rezultat je namreč eliprični cilinder, ki pa ima napačne koordinate centra.
Ali je to posledica tega, kako dobim matriko MM iz katere iščem člene G,H,I? Tam namreč upoštevam "star" center. Nevem kako drugače naj dobim člene G,H,I.
Image

Lp; Jan

User avatar
Aniviller
Posts: 7263
Joined: 15.11.2004 18:16

Re: Konstrukcija elipse

Post by Aniviller » 7.1.2010 13:31

Ce vse prav naredis, mora imeti popravljeni cilinder isti polosi in isti center. Ce ni tako, je nekje nekaj narobe.

Par pripomb:
Absolutna vrednost pri Q: Negativen Q pomeni, da ni resitev. Tvoja napaka je napacen preznak Q-ja, ker nisi uposteval, da smo nesli vse tiste konstante cez enacaj.
Ker so lambde sortirane, je lambda[[1]] in lambda[[2]] res prava izbira (zadnja bo najmanjsa po absolutni vrednosti). To je ok.

Code: Select all

P={{vector [[1, 1]], vector [[1, 2]], vector [[1, 3]]}, {vector [[2, 1]], vector [[2, 2]], vector [[2,3]]},
{vector [[3, 1]], vector [[3, 2]], vector [[3, 3]]}}
Hej! "vector" je ze P - zakaj prepisujes elemente enega po enega?
Za diagonalno matriko lahko uporabis DiagonalMatrix[{lambda[[1]],lambda[[2]],0}].
Elipticni cilinder ima seveda obe lastni vrednosti pozitivni (tretja je 0). Ce si pravo postavil na 0, mora izgledati enako. Vrstni red koeficientov mora biti isti - A,B,C,... so na istih mestih v matriki.

Zdaj pa pomemben del. Center. Ko si popravil matriko (M2 namesto M), mora biti center se vedno na istem mestu (x0,y0,z0 so isti). Spremenijo se pa G,H,I. Nove dobis nazaj iz iste enacbe iz katere si izracunal center.
Reseval si enacbo \(M\vec{r}_0=-(G/2,H/2,I/2)\). Nove torej dobis takole: \(-(G'/2,H'/2,I'/2)=M_2\vec{r}_0\) oziroma
\((G',H',I')=-2M_2\vec{r}_0\)

Zdaj imamo linearne clene - spremeni se pa tudi konstantni clen. Prej je bil 1 (oziroma v splosnem J). Zdaj je malo drugacen. Q je pa isti. Torej, prej je bilo
\(-Q=A x0^2+...+F y_0z_0+J\)
Zdaj je pa
\(-Q=A' x0^2+...+F' y_0z_0+J'\)
se pravi
\(J'=-Q-(A' x0^2+...+F' y_0z_0)\)

Mogoce ti bo malo pomagalo pojasnilo, kako do vsega tega pridemo, ce pa ni razumljivo pa tudi ni probem:
pojasnilo wrote:Vsi ti koraki so razumljivi, ce pogledas kaj smo naredili.
Imamo enacbo s koeficienti, ISCEMO pa obliko enacbe \((r-r_0)M(r-r_0)=Q\) (r=(x,y,z), r0=(r0,y0,z0)). Mnozenje matrike z vektorjem z obeh strani da skalar (po komponentah \(r_1 r_1 M_{11}+r_1 r_2 M_{12}+r_1 r_3 M_{13}+\cdots\), verjetno vidis vzorec)

Ce zmnozimo zgornjo enacbo, pride
\(\underbrace{rMr}_{\text{kvadratni cleni}}+\underbrace{2r Mr_0}_{\text{linearni cleni}}+\underbrace{r_0 M r_0-Q}_{J}=0\)

Torej, iz kvadratnih clenov smo prebrali M, iz linearnih clenov smo dobili tisto enacbo, iz katere smo dolocili center \(r_0\), zadnji clen nam je pa povedal Q.
(\(r_0 M r_0\) ni nic drugega kot tisti \(A x_0^2+By_0^2+\cdots\)).

M opisuje obliko in orientacijo ploskve, Q pa njeno velikost.

Ko smo popravili matriko, je samo pomenilo, da moramo novo matriko vstaviti v to enacbo namesto M, kar nam takoj da izraze, kako dobiti nove koeficiente.
Saj verjetno te bo v koncni fazi bolj zanimalo, kaksne so polosi in kaksna je smer in premik cilindra, tako da ne vem ce bo pridobivanje popravljenih koeficientov sploh potrebno, razen za risanje ploskve.

Medo5
Posts: 26
Joined: 8.2.2007 9:08
Location: Ljubljana
Contact:

Re: Konstrukcija elipse

Post by Medo5 » 7.1.2010 14:34

Aniviller wrote:Ce vse prav naredis, mora imeti popravljeni cilinder isti polosi in isti center. Ce ni tako, je nekje nekaj narobe.
Center res ostane enak. Ampak iz slike iz prejšenega odgovora je razvidno, da so moje točke ravno v nasprotju z narisanim cilindrom. Ali nebi bilo smiselno, da bi potekal cilinder skozi čim več točk na moji sliki? To bi bil najboljši fit. Sedaj pa, kot da ga ni.
Še nekaj sem opazil: če na zgorni sliki razpolovim z0 in zamenjam predznak torej z0=-z0/2 pride ravno popolen fit. Ampak mislim, da je to zgolj naključje.

Skratka zakaj je center izven mojh točk in ne proti sredini. Oz vem zakaj je, ker je bil prvi fit hiperbolični cilinder, ki pa ima center zunaj mojih točk(točno na tem mestu). Vprašanje je samo kako rešiti to.

Lp; Jan

User avatar
Aniviller
Posts: 7263
Joined: 15.11.2004 18:16

Re: Konstrukcija elipse

Post by Aniviller » 7.1.2010 14:53

Najprej poskusi vnest vse v zvezi s centrom, kar sem zapisal zgoraj in pogledat ce imas nove koeficiente v redu. Najbolj bi pa pomagalo, ce bi dal obe sliki: prvoten fit, in potem ko naredis cilindricnega (pa ce bi dal center na sliko kot se eno tocko). Potem se da hitro ugotovit ce je ze center narobe dolocen ali je kje drugje napaka.

Pa se nekaj: fit bo brezpogojno uspesen samo, ce bodo tocke bolje razporejene okrog cilindra. Ker na kratek lok, ki ga vidim na sliki, lahko nafitas skoraj karkoli - tudi hiperbolicen cilinder.

User avatar
Aniviller
Posts: 7263
Joined: 15.11.2004 18:16

Re: Konstrukcija elipse

Post by Aniviller » 7.1.2010 15:12

Aja, se tisto pri Q, ni cel predznak narobe... Q= Ax0^2+By0^2+...-j (oziroma -1 ce je pri tebi j=1). Tako bo pa ok. Sem preveril na hitro na Mathematici.

Medo5
Posts: 26
Joined: 8.2.2007 9:08
Location: Ljubljana
Contact:

Re: Konstrukcija elipse

Post by Medo5 » 7.1.2010 15:14

Tako, prilagam sliki:
na sliki je hiperbolični in eliptični cilinder. Oba imata isti center. Ralika ja ali posilim obe lambdi, da sta pozitivni. Sploh prvi fit je pa čisto ena tretja oblika. Seveda z istim centrom.(ni na sliki)
Image
Image
Pa se nekaj: fit bo brezpogojno uspesen samo, ce bodo tocke bolje razporejene okrog cilindra. Ker na kratek lok, ki ga vidim na sliki, lahko nafitas skoraj karkoli - tudi hiperbolicen cilinder.
Točke niso idealne, ravno zaradi tega bi jih rad popisal matematično. Enačbe pa "nastavil" tako, da pride ven eliptični cilinder.

Lp; Jan

User avatar
Aniviller
Posts: 7263
Joined: 15.11.2004 18:16

Re: Konstrukcija elipse

Post by Aniviller » 7.1.2010 15:36

Medo5 wrote:Točke niso idealne, ravno zaradi tega bi jih rad popisal matematično. Enačbe pa "nastavil" tako, da pride ven eliptični cilinder.
Saj ni treba, da so idealne, samo ce so preblizu skupaj ne povedo veliko o obliki.

Mozno je, da si narobe fital. Sem v Mathematici reproduciral cel postopek na kratek in pregleden nacin (razen fitanja, ce imas tam napako, potem bos moral sam ugotovit kaj je narobe), pa si poglej.
http://www.easy-share.com/1908937352/FitCil.nb

Medo5
Posts: 26
Joined: 8.2.2007 9:08
Location: Ljubljana
Contact:

Re: Konstrukcija elipse

Post by Medo5 » 7.1.2010 17:56

Sem poizkusil tvoje točke(enačbo) z mojim postopkom in rezultat je enak.
Moj problem->Ali se je možno že na začetku omejit, da iščem cilinder? Prilagam mojo datoteko iz Mathematice. Matematično je središče na pravem mestu, vendar drugje, kjer bi si želel.
Se da kaj postoriti?
http://www.easy-share.com/1908938286/Ci ... ticaJan.nb

Lp; Jan

User avatar
Aniviller
Posts: 7263
Joined: 15.11.2004 18:16

Re: Konstrukcija elipse

Post by Aniviller » 7.1.2010 18:22

Vidim ja... sicer ni cisto pravilno tisto na koncu (Q nimas prav popravljen in nista oba enako siroka). Problem je, da so tocke samo na eni strani cilindra in je resitev prakticno katerakoli kvadratna forma. Nekaj tock na drugi strani cilindra bi precej pomagalo.

Je mogoce omejit, da isces cilinder... vendar imas v tem primeru opravka z vezanim ekstremom, kar pomeni da vsa ta enostavnost z matrikami odpade... lahko pa uporabis vgrajene funkcije Mathematice. Saj hujsi problem kot cilindricnost je to, da mora biti elipticen (to ne bi bil problem, ce bi imel vsaj kaksno tocko na drugi strani). Bom malo pogledal, pa bom povedal ce se obnese.

User avatar
Aniviller
Posts: 7263
Joined: 15.11.2004 18:16

Re: Konstrukcija elipse

Post by Aniviller » 7.1.2010 19:28

Uh... fit z vezanim ekstremom totalno propade... resitve ne najde (pa se ni nobenega pametnega nacina za zagotovit pozitivne lastne vrednosti, kar je bistveno bolj pretec problem).
Lahko bi mogoce nasli se kak podoben nacin... ampak iskanje resitve bo vedno zelo nestabilno, ker se algoritem vedno ujame v lokalne resitve.

Zdaj ce bodo realni nesimulirani podatki taki, da bodo imeli se kaksno tocko na nasprotni strani, potem bo v redu. Tvoji trenutni podatki so pa zelo neprimerni za fitanje.

Sicer vidim se eno moznost, ki bi lahko delovala ampak bi jo bilo treba konkretno premislit. Namrec - lastne smeri dobimo prakticno pravilne (matrika P). Torej iscemo samo linearne clene in dve lastni vrednosti. Z drugimi besedami, namesto
\(Ax^2+By^2+Cz^2+Dxy+Exz+Fyz+Gx+Hy+Iz=1\)
imas potem enacbo
\(\lambda_1 (v_1\cdot (x,y,z))^2+\lambda_2 (v_2\cdot(x,y,z))^2+Gx+Hy+Iz=1\)
(notri nastopajo skalarni produkti s prvim in drugim lastnim vektorjem). Vendar tega tudi ne smemo resevat tako kot prej (tvoja matrika A, stolpec enk, psevdoinverz) ampak bi morali uporabiti FindFit[], pri cemer bi posebej povedali, da morata biti lambdi pozitivni.
FindFit bi imel za prvi argument podatke (pts, samo da ima vsaka komponenta se enko na koncu - za desno stran enacbe), drugi argument bi bil {l1 (v1* (x,y,z))^2+l2 (v2*(x,y,z))^2+Gx+Hy+Iz,lambda1>0,lambda2>0}, tretji {l1,l2,G,H,I} in cetrti seveda {x,y,z}. Nisem pa probal.

Medo5
Posts: 26
Joined: 8.2.2007 9:08
Location: Ljubljana
Contact:

Re: Konstrukcija elipse

Post by Medo5 » 7.1.2010 19:30

Točka na drugi strani cilindra res generira "tubo". Pač določim eno na pribljižno. Moti me le, ker ima strašen vpliv. Verjetno zato, ker je precej oddaljena od preostale gruče.

Imam pa problem z radiji. Ali si slučajno opazil, da je kaj narobe pri izračunu? Namesto 2.1 in 3.4 je rezultat 9.6 in 15.6.
Zanimivo pa je, da so vsi ostali parametri popolnoma pravilni!

Pri meni v Mathematici člen center.M.center ne deluje. Oz deluje samo na tvojih številih. Na mojih pa ne... Nevem v čem je problem. Napaka pa je:

Code: Select all

QQ = CENTER.M.CENTER
R1 = Sqrt[Abs[Q/lambda[[1]]]]
R2 = Sqrt[Abs[Q/lambda[[2]]]]
Dot::dotsh: Tensors {{5.34483},{5.34483},{8.10183*10^-16}} and \
{{-0.0130719,-0.00588235,-1.73472*10^-18},{-0.00588235,-<<21>>,-<<24>>\
},{-1.73472*10^-18,-1.13841*10^-18,-4.63496*10^-18}} have \
incompatible shapes. >>
Lp; Jan

User avatar
Aniviller
Posts: 7263
Joined: 15.11.2004 18:16

Re: Konstrukcija elipse

Post by Aniviller » 7.1.2010 20:00

Recimo QQ je ze narobe (tista plus enka ti manjka), kar bi lahko pripomoglo k napacnemu radiju. Dobro si poglej tisto kar sem ti poslal, da bos lahko s primerjavo nasel kaksne tezave.

Drugace je pa razlog za tole napako, da je pri tebi CENTER stolpec (v bistvu je matrika 3x1), ker si B definiral na ta nacin. B definiraj kot
B = {0.5, 1, 1.5};
(z enojnimi zavitimi oklepaji). Potem bo tudi center tako definiran in bo vse delovalo.

Tocka na drugi strani zato tako vpliva, ker imas brez nje definirano samo neposredno okolico temena na eni strani - od tam naprej pa lahko zavije kamor koli, ne da bi bistveno vplivalo na ujemanje v temenu (v tem primeru je zavilo navzven v hiperbolicno obliko). Ce dodas tocko na drugi strani, ga prisilis, da zaokrozi in zavije nazaj, da pobere to tocko. S to tocko takorekoc "zapnes" cilinder. Ce bos v podatkih imel tovrstne tocke, potem se lahko zaneses na enako stabilnost kot pri tisti 2d verziji kjer si fital elipso.

Medo5
Posts: 26
Joined: 8.2.2007 9:08
Location: Ljubljana
Contact:

Re: Konstrukcija elipse

Post by Medo5 » 8.1.2010 23:28

Radiji so bili napačni, ker sem pozabil en plus v izračunu Q-ja. Sedaj bo vredu. Vse skupaj že kodiram v C-ju. Ko bom prišel do končnega rezultata, ga bom objavil. :)

Lp; Jan

Post Reply