Hypertason sovittamisesta ja suurista käänteismatriiseista

Seuraa 
Viestejä1473
Liittynyt6.6.2009

Aluksi ajattelin kirjoittaa kutakuinkin tälläisen viestin tänne:

petsku
Oletetaan, että on kuutiollinen hila massa-alkioita, jotka ovat jousilla kytketyt naapureihinsa (unohdetaan kuitenkin viistot ja diagonaaliset naapurit), eli jokaisella on kuusi naapuria. Nämä alkiot voivat sitten värähdellä vapaasti naapuriensa jousivoimien valjastamina kolmessa ulottuvuudessa.
Jos on kolme pistettä, niin niihinhän voi kivasti sovittaa paraabelin ja interpoloida sen derivaatan avulla kulmakerrointa tällä välillä. Yleistin tuon saman periaatteen kolmeen ulottuvuuteen ja muodostin matriisisysteemin hypertason sovittamiseksi vastaavasti 3D-hilapisteisiin. Hypertasohan on periaatteessa tiheysfunktio 3D-avaruudessa, joten siitä pitäisi saada hyvähkö approksimaatio hiukkasten tiheysgradientille hilassa, eikö?
No matriisimuotoinen yhtälö tilanteestahan on:
Aw=k,
missä w pitää sisällään kaikkien kuuden naapurin - i, j, k, l, m, n - poikkeamat origoon asetetusta tarkastelualkiosta kaikkien kolmen akselin suhteen - x, y, z - siis:
w=(xi, yi, zi, xj, yj, zj, xk, yk, zk, xl, yl, zl, xm, ym, zm, xn, yn, zn)^T
ja k sulkee sisäänsä ax^2+by^2+cz^2+dx+ey+fz -muotoisen hypertason parametrit:
k=(a, b, c, d, e, f)^T
Tuo matriisi A onkin sitten astetta rumempi turilas:
[code:1599au32]/ 0 yi^2 zi^2 0 yi zi \
| xi^2 0 zi^2 xi 0 zi |
| xi^2 yi^2 0 xi yi 0 |
| 0 yj^2 zj^2 0 yj zj |
| xj^2 0 zj^2 xj 0 zj |
| xj^2 yj^2 0 xj yj 0 |
| 0 yk^2 zk^2 0 yk zk |
| xk^2 0 zk^2 xk 0 zk |
| xk^2 yk^2 0 xk yk 0 |
| 0 yl^2 zl^2 0 yl zl |
| xl^2 0 zl^2 xl 0 zl |
| xl^2 yl^2 0 xl yl 0 |
| 0 ym^2 zm^2 0 ym zm |
| xm^2 0 zm^2 xm 0 zm |
| xm^2 ym^2 0 xm ym 0 |
| 0 yn^2 zn^2 0 yn zn |
| xn^2 0 zn^2 xn 0 zn |
\ xn^2 yn^2 0 xn yn 0 /[/code:1599au32]
Täällä on varmasti joku, joka on joskus jotakin vastaavaa puuhastellut ja osaisi sanoa, onko tuo matriisini aivan hakuteillä. Yleistys paraabelistä tason kautta hypertasoon on kuitenkin varsin suoraviivainen, joten toivottavasti se on edes hiukan sinnepäin.
Tämän matriisiyhtälön saa mukavasti ratkaistua PNS:llä:
A^T Aw=A^T k
Tai mukavasti ja mukavasti - jo A^T A -matriisitulon muodostaminen on vähintäänkin paskamaista käsipelillä, mutta sitä varten voin koodata pienen ohjelman. Lopulta minun kuitenkin pitäisi ratkaista 18x18-käänteismatriisi. On sanomattakin selvää, että ulkoistan sen tietokoneelle, mutta mitä matemaattista vippaskonstia minun kannattaisi käyttää? Ehkä LU-hajotelma olisi oikotie onneen? Osaisiko joku vihjata jotakin inhimillistä keinoa?

Mutta totesin sitten tätä kirjoittaessani, että matriisini kuvaa korkeintaan poikittaisen aaltoliikekomponentin tiheyttä, jos sitäkään. Osaisiko joku, jolla on matemaattista silmää, sanoa, mitä ihmettä se oikein kuvaa? Tajusin myös jo, että kaikki, mitä oikeasti tarvitsen, on pitkittäissuuntaisten poikkeamien tarkastelu ja miten se toteutetaan miljoona kertaa yksinkertaisemmin.
No tuli tosiaan pitkä viesti kirjoitettua, niin ei sitä raaski poistaakaan - yleistä keskustelua em. ja niitä sivuavista aiheista, mikäli ne ketään muuta kiinnostavat, kiitos, hei.

Kommentit (7)

Neutroni
Seuraa 
Viestejä26890
Liittynyt16.3.2005
petsku
. Lopulta minun kuitenkin pitäisi ratkaista 18x18-käänteismatriisi. On sanomattakin selvää, että ulkoistan sen tietokoneelle, mutta mitä matemaattista vippaskonstia minun kannattaisi käyttää? Ehkä LU-hajotelma olisi oikotie onneen? Osaisiko joku vihjata jotakin inhimillistä keinoa?



Tarvitsetko käänteismatriisin vai vain ratkaisun yhtälöryhmälle? Onko sinulla jokin kiire? 18 yhtälön lineaarinen ryhmä ratkeaa ihan tavallisella Gaussin algoritmilla ja kääntyy Gaussilla ja Jordanilla. Mutta älä käännä sitä, jos et tarvitse sitä montaa kertaa.

Osaisiko joku, jolla on matemaattista silmää, sanoa, mitä ihmettä se oikein kuvaa?



Enpä jaksa alkaa miettiä. Miksi asiaa ylipäätään pitää lähestyä noin monimutkaisesti? Lasket vaan ensin hiukkasiin kohdistuvat kiihtyvyydet ja niistä edelleen uudet nopeudet ja paikat. Toistat sitä tarpeeksi, niin saat toivomasi tuloksen.

No tuli tosiaan pitkä viesti kirjoitettua, niin ei sitä raaski poistaakaan - yleistä keskustelua em. ja niitä sivuavista aiheista, mikäli ne ketään muuta kiinnostavat, kiitos, hei.



Onko tuosta jotain hyötyä vai huviksesiko simuloit virtuaalista aladoobia?

petsku
Seuraa 
Viestejä1473
Liittynyt6.6.2009
Neutroni
Tarvitsetko käänteismatriisin vai vain ratkaisun yhtälöryhmälle? Onko sinulla jokin kiire? 18 yhtälön lineaarinen ryhmä ratkeaa ihan tavallisella Gaussin algoritmilla ja kääntyy Gaussilla ja Jordanilla. Mutta älä käännä sitä, jos et tarvitse sitä montaa kertaa.

Eipä taas ole hetkeen ollut matriisimatikkaa vastassa, kun perus-gaussaaminen unohtui... Kiitos muistutuksesta!
Niin ja minkään valtakunnan kiirettä ei ole.

Neutroni
Enpä jaksa alkaa miettiä. Miksi asiaa ylipäätään pitää lähestyä noin monimutkaisesti? Lasket vaan ensin hiukkasiin kohdistuvat kiihtyvyydet ja niistä edelleen uudet nopeudet ja paikat. Toistat sitä tarpeeksi, niin saat toivomasi tuloksen.

Siis itse hiukkasten simuloiminen on kyllä hanskassa ja sen toteutan RK4:llä, mutta olisin halunnut kehittää työkalun visualisointiin, ja kuten itsekin totesin, niin metsään meni ja syvälle.
Lisäksi olisin saanut tekosyyn koodata itselleni paketit vektori- ja matriisioperaatioita varten.
Neutroni
Onko tuosta jotain hyötyä vai huviksesiko simuloit virtuaalista aladoobia?

Ihan huvikseni vain - huvinsa kullakin.
Nyt kun on taas aikaa, niin voisi tehdä ohjelman värähtelylle kaikissa kolmessa ulottuvuudessa ja tutkia vaikka akustiikkaa ja mahdollisesti törmäyksiä, jolloin kimmokertoimen voisi korvata jännitys-venymä-käyrällä murtorajoineen. Plastisuudenkin uskoisin osaavani alkeellisesti ympätä tuohon. Nämä nyt voisi toteuttaa vaikka COMSOLilla, jonka parissa tullee nyt kesällä rutkasti aikaa vietettyä, mutta matematiikka näiden takana ja sen keksiminen itse kiinnostaa ja harva asia on mielenkiintoisempaa, kuin leikkiä jumalaa ja luoda omat maailmansa.

Vierailija
petsku
Oletetaan, että on kuutiollinen hila massa-alkioita, jotka ovat jousilla kytketyt naapureihinsa (unohdetaan kuitenkin viistot ja diagonaaliset naapurit), eli jokaisella on kuusi naapuria. Nämä alkiot voivat sitten värähdellä vapaasti naapuriensa jousivoimien valjastamina kolmessa ulottuvuudessa.
Jos on kolme pistettä, niin niihinhän voi kivasti sovittaa paraabelin ja interpoloida sen derivaatan avulla kulmakerrointa tällä välillä. Yleistin tuon saman periaatteen kolmeen ulottuvuuteen ja muodostin matriisisysteemin hypertason sovittamiseksi vastaavasti 3D-hilapisteisiin. Hypertasohan on periaatteessa tiheysfunktio 3D-avaruudessa, joten siitä pitäisi saada hyvähkö approksimaatio hiukkasten tiheysgradientille hilassa, eikö?

Jäi nyt kyllä epäselväksi, miksi haluat sovittaa paraabelia kolmeen avaruuden pisteeseen ja ovatko nuo kolme pistettä hiukkasten koordinaatteja? Voisitko vähän valaista, mitä tuo paraabeli taikka sen kulmakerroin on kuvaavinaan?

Jos katsotaan hiukkastiheyden määritelmää, niin sehän on yhdessä ulottuvuudessa dn(x)/dx ja kolmessa ulottuvuudessa dn(x)/dV, siis hiukkaslukumäärän n(x) derivaatta tilavuuden suhteen. Tiheyden gradientti taas on d²n(x)/dV².

Jos siis haluat laskea hiukkasssysteemisi hiukkastiheysjakauman (tai massatiheysjakauman) niin sinun pitää ensin jakaa alueesi sopiviin koppeihin ja katsoa mikä on hiukkaslukumäärä n(x). Tiheys selviää sitten derivoimalla sitä tilavuuden suhteen ja gradientti tupladerivoimalla.

Vanha jäärä
Seuraa 
Viestejä1557
Liittynyt12.4.2005
petsku

Tämän matriisiyhtälön saa mukavasti ratkaistua PNS:llä:
A^T Aw=A^T k

Olen joskus askarrellut tällaisten sovitusten kanssa - en tosin mitenkään suurilla muuttujamäärillä, ja siksi varoittaisin tällaisesta suorasta lähestymistavasta. Vähänkään korkea-asteisempia polynomeja sovitettaessa on vaarana merkitsevien numeroiden ylivuoto, mikä tekee ratkaisusta mitä sattuu. Näin tapahtuu ainakin silloin, kun sovitettava parametrialue on laaja.

Sitä vastoin olisin ratkaisemassa perusyhtälöä Aw=k SV-hajotelman avulla. Tiedän kyllä, että SVD on numeerisesti raskas, mutta se kokemukseni mukaan tuottaa aina tarkat ja stabiilit ratkaisut.

Vanha jäärä

petsku
Seuraa 
Viestejä1473
Liittynyt6.6.2009
Schlierenzauer
Jäi nyt kyllä epäselväksi, miksi haluat sovittaa paraabelia kolmeen avaruuden pisteeseen ja ovatko nuo kolme pistettä hiukkasten koordinaatteja? Voisitko vähän valaista, mitä tuo paraabeli taikka sen kulmakerroin on kuvaavinaan?

Siis se, mitä hain takaa ja minkä olinkin ensin vähän nurinkurin ajatellut, on kutakuinkin seuraavaa:
Tarkastellaan kolmea hiukkasta x-akselin suhteen, jotka ovat toistensa suhteen levossa pisteissä -1, 0 ja 1. Kun hiukkasta poikkeuttaa näistä lepoasemistaan, eli niiden välinen etäisyys ≠ 1, kokevat ne jousivoiman. Valitaan tuo koordinaatisto siten, että keskimmäinen hiukkanen pysyy aina origossa. Sanotaan nyt vaikka, että reunimmaiset hiukkaset päätyvät keskimmäisen suhteen paikkoihin -0,8 ja 1,3. Tällöin niiden venymät ovat -0,8 ja 1,3. Sijoitetaan venymät venymä-akselille e, joka on x-akseliin nähden kohtisuorassa. Käytetään nyt x-koordinaatteina lepoasemia ja e-koordinaatteina venymän arvoja. Saadaan pisteet (-1; -0,8), (0, 0) ja (1; 1,3). Kun näihin pisteisiin sovitetaan paraabeli, saadaan käsittääkseni approksimaatio tiheysgradientin x-komponentille tuon keskimmäisen hiukkasen kohdalla.
Minä vain aluksi lähdin jostakin ihmeen syystä pyörittelemään poikittaisten siirtymien kanssa ja eihän se putkeen mennyt.

E: Tarkemmin ajatellen se taitaakin olla enemmänkin harventumagradientti... Osoittaa siis harvemman tavaran suuntaan... ehkä.

Neutroni
Seuraa 
Viestejä26890
Liittynyt16.3.2005

Jos haluat tiheydelle arvon jokaisen pisteen kohdalla, voit laskea tuolla tavalla suorakulmaisen särmiön vetämällä tasot aina naapurien puoliväliin. Särmiön tilavuus on kääntäen verrannollinen tiheyteen. Tuo toimii, jos pisteiden poikkeamat tasapainoasemasta ovat pieniä.

myl
Seuraa 
Viestejä224
Liittynyt18.11.2010
petsku
Neutroni
Neutroni
Nyt kun on taas aikaa, niin voisi tehdä ohjelman värähtelylle kaikissa kolmessa ulottuvuudessa ja tutkia vaikka akustiikkaa:twisted:



Tässä siis yritetään kehittää elastisuuden teoriaa.
Kannattaa varmaan luntata hieman vinkkejä alan kirjallisuudesta.
Aiheen kehittely on sentään vienyt huippufyysikoilta pari sataa vuotta.

Lyhyesti
jännitys = jäykkyys * venymä
ja nämä kaikki ovat tensoreita, jotka tosin muunnetaan matriiseiksi ja vektoreiksi (Voigtin muunnoksella) ennen laskemista.
Värähtelyongelmat johtavat Helmholtzin yhtälöön, jonka ratkaisu vaatii Christoffelin matriisin ominaisarvojen määrärämistä.

-myl

Uusimmat

Suosituimmat