E¨ otv¨ os Lor´ and Tudom´ anyegyetem Informatikai Kar Numerikus Anal´ızis Tansz´ ek
H´ uros hangszer modellje ´ es annak implement´ al´ asa
Stoyan Gisbert
Potempski D´ aniel Konr´ ad
Egyetemi tan´ ar
Programtervez˝ o Informatikus MSc
Budapest, 2012
Tartalomjegyz´ ek 1. Bevezet´ es
3
1.1. Hangszerek fizikai modellez´es´enek n´eh´any alkalmaz´asa . . . . . . . . .
4
1.2. Technikai megjegyz´es . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
2. A zongora v´ azlatos fel´ ep´ıt´ ese, m˝ uk¨ od´ ese
6
3. Az egyes o ¨sszetev˝ ok modelljei
8
3.1. A h´ ur modellj´evel kapcsolatos eddigi n´eh´any eredm´eny . . . . . . . .
8
3.2. A rezon´ans modellj´enek fel´ır´asa . . . . . . . . . . . . . . . . . . . . . 11 3.2.1.
A rezon´ans feladata ´es fel´ep´ıt´ese . . . . . . . . . . . . . . . . 11
3.2.2.
Eddigi eredm´enyek a fizikai modellek kutat´as´aban . . . . . . 11
3.2.3.
Az egydimenzi´os modell . . . . . . . . . . . . . . . . . . . . . 12
3.2.4.
A k´etdimenzi´os modell
. . . . . . . . . . . . . . . . . . . . . 13
4. Az o ¨sszetev˝ ok o ¨sszekapcsol´ as´ anak modelljei
15
4.1. A h´ ur ´es a rezon´ans k¨ozvetlen ¨osszekapcsol´asa . . . . . . . . . . . . . 15 4.1.1. Az egydimenzi´os eset . . . . . . . . . . . . . . . . . . . . . . . 15 4.1.2. A k´etdimenzi´os eset . . . . . . . . . . . . . . . . . . . . . . . . 15 4.2. A h´ ur ´es a rezon´ans o¨sszekapcsol´asa egy egyszer˝ u h´ıdon kereszt¨ ul . . 16 4.3. Az egyszer˝ u h´ıdon kereszt¨ ul val´o ¨osszekapcsol´as egy v´altozata . . . . 17 4.4. A gerjeszt´es modellez´ese . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.4.1. Kezdeti ´ert´ekkel val´o gerjeszt´es . . . . . . . . . . . . . . . . . 17 4.4.2. A jobboldalon kereszt¨ ul t¨ort´en˝o gerjeszt´es . . . . . . . . . . . 18 5. Numerikus megval´ os´ıt´ asok
19
5.1. A h´ ur egy egyszer˝ u modellj´enek megval´os´ıt´asa . . . . . . . . . . . . . 19 5.1.1. A v´egesdifferencia m´odszer fel´ır´asa . . . . . . . . . . . . . . . 19 5.1.2. Stabilit´as . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 5.1.3. Tesztek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 5.2. Az egydimenzi´os rezon´ans egyszer˝ u modellj´enek megval´os´ıt´asa . . . . 22
5.2.1. A v´egesdifferencia m´odszer fel´ır´asa . . . . . . . . . . . . . . . 22 5.2.2. Stabilit´as . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 5.2.3. A m´odszer tesztel´ese . . . . . . . . . . . . . . . . . . . . . . . 24 5.3. A k´etdimenzi´os rezon´ans egyszer˝ u modellj´enek megval´os´ıt´asa . . . . . 24 5.3.1. A v´egesdifferencia s´ema fel´ır´asa . . . . . . . . . . . . . . . . . 24 5.3.2. Stabilit´as . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 5.3.3. A m´odszer tesztel´ese . . . . . . . . . . . . . . . . . . . . . . . 27 5.4. A biharmonikus oper´ator saj´at´ert´ek probl´em´aj´anak numerikus megold´asa befogott v´ekony izotr´op lemez peremfelt´etelei mellett . . . . . 31 5.5. A k¨ozvetlen o¨sszekapcsol´as megval´os´ıt´asa . . . . . . . . . . . . . . . . 32 5.5.1. Az egydimenzi´os eset . . . . . . . . . . . . . . . . . . . . . . . 32 5.5.2. Az egydimenzi´os eset tesztel´ese . . . . . . . . . . . . . . . . . 34 5.5.3. A k´etdimenzi´os eset . . . . . . . . . . . . . . . . . . . . . . . . 36 5.6. A k´etdimenzi´os eset tesztel´ese . . . . . . . . . . . . . . . . . . . . . . 36 5.7. Az egyszer˝ u h´ıdon kereszt¨ uli o¨sszekapcsol´as megval´os´ıt´asa . . . . . . . 38 5.7.1. A s´ema ismertet´ese . . . . . . . . . . . . . . . . . . . . . . . . 38 5.7.2. Tesztel´es . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5.8. Az egyszer˝ u h´ıdon kereszt¨ uli ¨osszekapcsol´as megval´os´ıt´as´anak v´altozata 39 5.8.1. A s´ema ismertet´ese . . . . . . . . . . . . . . . . . . . . . . . . 39 5.8.2. Tesztel´es . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 6. A tov´ abbi kutat´ omunka lehets´ eges ir´ anyai
42
6.1. Pontosabb fizikai modellek haszn´alata . . . . . . . . . . . . . . . . . . 42 6.2. Pontosabb numerikus modellek, tesztek haszn´alata
. . . . . . . . . . 42
7. Felhaszn´ alt technol´ ogi´ ak, er˝ oforr´ asok, ´ es az azokkal kapcsolatos tapasztalatok 8. Eredm´ enyek o ¨sszefoglal´ asa
44 46
1. fejezet Bevezet´ es Jelen munka bizonyos hangszerek egyszer˝ u fizikai modellj´et ´ırja fel, majd ismerteti a numerikus megval´os´ıt´as mik´entj´et. A fizikai modelleket ismertetem, r¨oviden vitatom a benn¨ uk rejl˝o lehet˝os´egeket. K¨ozl¨ok n´eh´any lehets´eges sz´am´ıt´og´epes megval´os´ıt´ast, azok hib´ait, eredm´enyeit, ´es ezeket is vitatom. A tov´abbfejleszt´es t¨obb lehet˝os´eg´evel z´arok. A kidolgoz´as sor´an els˝odlegesen a zongor´ara f´okusz´altam, de a modell a´tvihet˝o a legt¨obb h´ uros hangszerre. A megval´os´ıt´asok, tesztek sor´an keletkezett programok megtal´alhat´ok a mell´ekelt adathordoz´on. A Diplomamunka-t´ema bejelent˝o lapon kit˝ uz¨ott t¨obb c´elt nem val´os´ıtottam meg, m´asokat viszont r´eszletesen, t¨obb megold´asi lehet˝os´eget is bemutatva dolgoztam ki. A hangszer modelljei fel´all´ıt´asakor nem vettem figyelembe a leveg˝ot ´es a nemline´aris effektusokat, ugyanis a szakirodalmat kutatva kider¨ ult, hogy elfogadott mind k´et jelens´eg elhanyagol´asa. Ez k¨ ul¨on¨osen igaz a h´ urok egym´asra hat´as´ara a leveg˝on kereszt¨ ul. A nemline´aris jelens´egekre viszont igaz, hogy van l´etjogosults´aguk [1], de az ennyire pontos modellek meghaladj´ak eme diplomamunka kereteit. A mell´ekelt programban a kett˝on´el t¨obb h´ ur eset´et is implement´altam, ami azt mutatja, hogy az ismertetett modellek term´eszetszer˝ uleg terjeszthet˝ok ki erre az esetre. A modellek lehets´eges diszkretiz´aci´oira ´es azok tulajdons´agaira ford´ıtottam nagyobb hangs´ ulyt. Nem foglalkoztam a modellparam´eterek meghat´aroz´as´aval, mert ez legink´abb a leveg˝o feljebb eml´ıtett hat´asakor lett volna relev´ans. M´er´esi eredm´enyekkel val´o o¨sszevet´es nem t¨ort´ent, viszont analitikus eredm´enyekkel igen. Csak a vizu´alis megjelen´ıt´es probl´em´aj´at oldottam meg. Mindenek el˝ott azonban felmer¨ ul a k´erd´es, hogy van-e gyakorlati haszna az ilyen modelleknek.
3
1.1.
Hangszerek fizikai modellez´ es´ enek n´ eh´ any alkalmaz´ asa
A v´alasz term´eszetesen igen. T¨obb ter¨ uleten is alkalmazhat´ok a modellek numerikus megval´os´ıt´asukkal egy¨ utt. El˝osz¨or is a hangszerek m˝ uk¨od´es´enek megismer´ese nagyban seg´ıtheti fejleszt´es¨ uket. A legt¨obb hangszer t¨obb sz´az ´eves fejl˝od´es nyom´an ´erte el jelenlegi form´aj´at, legink´abb sok-sok pr´ob´algat´as u ´tj´an fejlesztett´ek o˝ket. Igaz, volt amikor egyes r´eszletek meg´allap´ıt´as´ahoz (mint p´eld´aul f´ uv´os hangszerek lyukainak elhelyezked´ese) komoly sz´am´ıt´asokat v´egeztek. Az ilyen sz´am´ıt´asok elterjed´ese sokat seg´ıt a gy´art´oknak. Nishiguchi szerint addig nem ´all´ıthatjuk, hogy tiszt´an ´ertj¨ uk a zongora hangj´anak keletkez´es´et, am´ıg nem szimul´altuk azt hihet˝oen [27]. Ha ezt elfogadjuk, akkor a numerikus megval´os´ıt´asnak t¨obbsz¨or¨os szerepe van a hangszerek fejleszt´es´eben, ´ıgy az a´ltalunk vizsg´alt h´ uros hangszerek eset´eben is ez a helyzet. Nem csak az elm´eleti modell helyess´eg´et lehet ´ıgy bizony´ıtani, de k¨onnyebben kimutathat´o az egyes fizikai param´eterek megv´altoztat´as´anak hat´asa is a hangra. Kifejezetten a gy´art´ast seg´ıt˝o eredm´enyeket mutat be p´eld´aul A. Stulov [28]-ban. K´ıs´erletei sor´an numerikus szimul´aci´okat is alkalmazott. A Tallini Zongoragy´ar a javasolt v´altoztat´asokat meg is val´os´ıtotta. M´asik fontos alkalmaz´asi ter¨ ulet a fizikai modell alap´ u hangszint´ezis. A manaps´ag elterjedt jel alap´ u szint´ezissel szemben nagy el˝onye [1] alapj´an, hogy a modellparam´eterek az eredeti hangszer tulajdons´agait t¨ ukr¨ozik. ´Igy a felhaszn´al´o sz´am´ara sokkal k´ezzelfoghat´obb, hogy min is v´altoztat, ha m´odos´ıtja a param´etereket. R´aad´asul ezekb˝ol a param´eterekb˝ol kevesebb is el´eg, mert az alkalmazott modell bizonyos korl´atokat szab. H´atr´anya az ilyen megk¨ozel´ıt´esnek a nagy sz´am´ıt´asig´eny, ami egy val´os idej˝ u rendszer eset´en kompromisszumokhoz vezet. Sajnos a val´osidej˝ us´eg megtart´asa a modell egyszer˝ us´ıt´es´ehez, elnagyol´as´ahoz, s ´ıgy a gener´alt hang az eredetit˝ol val´o t´avolod´as´ahoz vezet. Lehets´eges azonban nem l´etez˝o hangszerek virtu´alis ´ep´ıt´ese is, ´es ezeket haszn´alni u ´j hangok szintetiz´al´as´ahoz. Egy fizikai modellen alapul´o mod´alis szint´ezist megval´os´ıt´o keretrendszer a MOSAIC [25]. Ennek alap¨otlete, hogy minden rezg˝o szerkezet ¨osszekapcsolhat´o mod´alis adatokkal. Ez persze csak akkor igaz, ha a modellek mind line´arisak. Teh´at minden rezg˝o szerkezet egy mod´alis objektum. Ezeket az objektumokat o¨ssze lehet k¨otni t¨obbf´ele m´odon is. A rendszerhez tartoznak m´eg gerjeszt˝o objektumok valamint virtu´alis pick-upok. Egy m´asik m´odszere az u ´j hangsz´ınek el˝o´all´ıt´as´anak [20]-ban van ismertetve. Ennek l´enyege, hogy k´et, fizikai modellen alapul´o hangsz´ın k¨oz¨otti interpol´aci´ot 4
v´egeznek. Az interpol´aci´o k´et m´odon mehet v´egbe. Az egyik eset, amikor a k´et hangsz´ın ugyanazzal a modellel van el˝oa´ll´ıtva, csak a modellparam´eterek k¨ ul¨onb¨oznek. Ekkor egyszer˝ uen a sz´oban forg´o param´etereket kell interpol´alni. A m´asik eset, amikor a modellek is k¨ ul¨onb¨oznek. Ekkor kell tal´alni egy harmadik, ´altal´anosabb modellt, amely mindkett˝onek a kiterjeszt´ese. Majd ezen bel¨ ul a modellen bel¨ ul kell megtal´alni az eredeti modelleknek megfelel˝o param´etereket, ´es azokat interpol´alni. A tanulm´any egy gyakorlati p´eld´an mutatja be a m´odszert. A git´ar ´es a zongora hangja k¨oz¨ott hoznak l´etre a szerz˝ok egy egyenletes a´tmenetet. A k´et hangszer modellje csak a gerjeszt´es m´odj´aban k¨ ul¨onb¨ozik, ez´ert csak azt kellett a´ltal´anos´ıtani. A m´odszer ´erdekess´ege, hogy az extrapol´aci´o is lehets´eges.
1.2.
Technikai megjegyz´ es
Mivel az a´ltalam hivatkozott szakirodalom nagy r´esze angol nyelv˝ u, ez´ert a fontosabb szakszavak els˝o felt˝ un´esekor z´ar´ojelben k¨ozl¨om az angol nyelv˝ u megfelel˝ot.
5
2. fejezet A zongora v´ azlatos fel´ ep´ıt´ ese, m˝ uk¨ od´ ese A zongoraj´at´ekos a billenty˝ uk lenyom´as´aval egy bonyolult mechanika k¨ozvet´ıt´es´evel hozza mozg´asba a kalap´acsot. A kalap´acs u ¨ti meg a h´ urt. A h´ urt a h´ıd (bridge) k¨oti ¨ossze a rezon´anssal (soundboard). A rezon´ans a zongora kem´enyfa test´ebe van be´agyazva. A h´ urok egyik v´eg´et a t˝ok´ebe (pin block) be´agyazott hangol´osz¨ogek (pin), m´asik v´eg´et a p´anc´elt˝oke (cast-iron frame) kiakaszt´o sz¨ogei (hitch pin) fesz´ıtik. A h´ ur egy ac´elmagb´ol (core) ´all. A m´elyebb h´ urok egy-, vagy k´etr´eteg˝ u r´ezfonatot (copper wrapping) kapnak. Az angol nyelv˝ u elnevez´eseket [12]-b˝ol vettem.
6
2.1. a´bra. A modern zongora f˝obb alkatr´eszei: 1-p´anc´elt˝oke, 2-fed´el(el¨ uls˝o), 3tomp´ıt´o, 4-tomp´ıt´o, 5-fed´el, 6-tomp´ıt´okar, 7-sostenuto s´ın, 8,9-ped´almechanizmus, 10-ped´alr´ ud, 11-ped´alok, 12-h´ıd, 13-kiakaszt´o sz¨og, 14-p´anc´elt˝oke, 15-rezon´ans, 16h´ ur. Forr´as: http://www.wikipedia.hu
7
3. fejezet Az egyes o ¨sszetev˝ ok modelljei ´ Ebben a fejezetben bemutatom a hangszer f˝o ¨osszetev˝oinek modelljeit. Altal´ anoss´agban h´arom t´ıpus´ u elemb˝ol a´ll egy h´ uros hangszer. H´ urokb´ol, a h´ urok vibr´al´as´at er˝os´ıt˝o eszk¨ozb˝ol, ami legt¨obbsz¨or falemez (zongora, csembal´o) vagy f´ab´ol k´esz¨ ult rezon´al´o doboz (git´ar, heged˝ u, h´arfa), de lehet membr´an is (banjo). Mivel a zongora modellez´ese az els˝odleges c´elunk, ez´ert a falemez eset´et vizsg´aljuk. Sz¨ uks´eges m´eg a h´ urt gerjeszteni valamilyen m´odon. Ez a harmadik elem, amely n´eha maga a j´at´ekos. Mi mag´at a gerjeszt´est fogjuk modellezni, nem a gerjeszt˝o eszk¨ozt.
3.1.
A h´ ur modellj´ evel kapcsolatos eddigi n´ eh´ any eredm´ eny
A h´ ur modellje manaps´ag a legapr´obb r´eszletekig kidolgozott. A legf˝obb eredm´enyek ´ [1, 12, 26] alkot´asokat veszem most alat¨obb m˝ uben is ¨ossze vannak foglalva. En pul. A h´ ur legegyszer˝ ubb megk¨ozel´ıt´ese r´eg´ota ismert. Ez egy idealiz´alt eset, ide´alis h´ urnak is szokt´ak nevezni. Az ilyen h´ ur mozg´as´at nem tomp´ıtja semmi, nem veszik el az energi´aja, nincs merevs´ege, azaz t¨ok´eletesen rugalmas, valamint csak az egyik mer˝oleges ir´anyba val´o kit´er´est vizsg´aljuk. Ebben az esetben a h´ ur mozg´asegyenlete a j´ol ismert hull´amegyenlet [12]: 2 ∂ 2u 2∂ u = c dt2 ∂x2
ahol c =
q
Θ , µ
Θ a fesz´ıt´es m´ert´eke, µ a t¨omeg osztva az egys´egnyi hosszal (vagy
t¨omegs˝ ur˝ us´eg), u(x, t) : [0, L] × [0, T ] → R (0 < L, T ∈ R) a keresett megold´as, azaz a kit´er´es a hely ´es az id˝o f¨ uggv´eny´eben. Levezet´ese megtal´alhat´o t¨obbek k¨oz¨ott ´ ett˝ol elt´erek, hogy [1, 12, 26]-ban. A szakirodalom Θ-t egys´egesen T -vel jel¨oli. En az id˝ointervallum v´eg´et jel¨olhessem vele. Ha a h´ urra hat´o k¨ uls˝o er˝ot is ki szeretn´enk 8
fejezni, akkor az inhomog´en µ
∂ 2u ∂ 2u − Θ =f dt2 ∂x2
(3.1)
egyenletre jutunk, ahol f (x, t) : [0, L] × [0, T ] → R, adott [1]. Peremfelt´etelk´ent adottak: u(0, t), u(L, t) (t ∈ [0, T ]). A kezdeti ´ert´ekek pedig: u(x, 0),
∂u (x, 0) (x ∈ [0, L]). ∂t
Ha a peremfelt´etelek lek¨ot¨ott v´egeknek felelnek meg, azaz u(0, t) ≡ 0,
u(L, t) ≡ 0,
(3.2)
akkor a saj´atfrekvenci´ak a k¨ovetkez˝ok[12, 1]: ωk = k
cπ L
(k ∈ N+ )
sz¨ogfrekvenci´ak eset´en, ´es c (k ∈ N+ ) 2L a kit´er´esek frekvenci´aja eset´en. Amikor a felhangok frekvenci´ai ilyen eg´esz sz´am´ u fk = k
t¨obbsz¨or¨os¨okb˝ol ´all´o sorozatot alkotnak, akkor harmonikusaknak h´ıvjuk o˝ket. Ez a tulajdons´ag a zen´eben fontos, hiszen a hangzatok (akkordok), ´es ´ıgy a kialakult sk´al´ak alapj´aul szolg´al, emiatt a hangszerek hangol´as´an´al is szerepe van. Viszont amint k¨ozel´ıtj¨ uk a h´ ur modellj´et a val´os´aghoz, a saj´atfrekvenci´ak kisebb-nagyobb m´ert´ekben el fognak t´erni a harmonikusokt´ol. Ezt h´ıvja a szakirodalom inharmonicit´asnak. A legkor´abban felfedezett tulajdons´ag, ami inharmonicit´ast okoz, a val´os h´ ur merevs´ege (stiffness). Az ilyen h´ ur kicsit r´ udk´ent is viselkedik, azaz nem csak a fesz´ıt´es´eb˝ol, hanem a merevs´eg´eb˝ol is ad´odik egy, az egyens´ ulyi a´llapotba visszat´er´ıt˝o er˝o. Az ilyen tulajdons´agot kifejez˝o negyedrend˝ u deriv´alt szerepel a r´ ud vagy ge2
renda (3.3) (idealiz´alt) mozg´asegyenlet´eben, (ahol majd a −Θ ∂∂xu2 fesz´ıt´est kifejez˝o m´asodrend˝ u tag nem is szerepel). Ezt a negyedrend˝ u deriv´altat alkalmazva a µ
4 ∂ 2u ∂ 2u 2∂ u − Θ + ESκ =f dt2 ∂x2 ∂x4
mozg´asegyenletre jutunk, ahol µ = ρS, ρ a s˝ ur˝ us´eg, S a keresztmetszet ter¨ ulete, E a Young-f´ele modulus, κ a forg´as sugara (radius of gyration [26]). A peremfelt´etelek is v´altoznak: a peremen is sz¨ uks´eges megadni deriv´altat a gerenda eset´ehez hasonl´oan. A f¨ uggesztett esetnek megfelel˝o peremfelt´etelek mellett 9
(p´aros helyderiv´altak null´ak, a nulladikat is bele´ertve) a k¨ovetkez˝o saj´atfrekvenci´ak ad´odnak [12, 1]:
√ fk = kf0 1 + Bk 2
(k ∈ N+ )
ahol B = π 2 ESκ2 /ΘL2 . Befogott v´egek eset´ere ((3.2) ´es az els˝orend˝ u deriv´altak is null´ak) [12, 26] k¨oz¨olnek approxim´aci´okat, melyek azt mutatj´ak, hogy a fenti k´eplethez k´epest l´enyeg´eben egy konstans szorz´o a k¨ ul¨onbs´eg. A merevs´eg nem felt´etlen homog´en, hisz a zongorah´ urok eset´eben a r´eztekercs bizonyos h´ urokn´al k´etr´eteg˝ u, ´es a bels˝o kor´abban kezd˝odik, mint a k¨ uls˝o. Teh´at a h´ ur ekkor h´arom k¨ ul¨onb¨oz˝o merevs´eggel ´es t¨omegs˝ ur˝ us´eggel rendelkezhet. Az inharmonicit´as egyr´eszt megnehez´ıti a zongora hangol´as´at. Az egyes hangk¨oz¨ok a magasabb hangok fel´e haladva egyre t´avolodnak egym´ast´ol. M´asr´eszt az inharmonicit´as d´ us´ıtja a spektrumot, egyesek ezt u ´gy fogalmazz´ak meg, hogy ,,melegebb´e” teszi a hangz´ast, azaz jav´ıtja a hangsz´ınt. Ennek hallhat´os´ag´at vizsg´alja [21], ´es azt ´allap´ıtja meg, hogy a m´elyebb hangok eset´en az ember sz´am´ara k¨onnyebben ´eszlelhet˝o a jelens´eg. A h´ urok s´ uly´anak n¨ovel´es´et a m´elyebb hangokn´al pont az inharmonicit´as hat´as´anak minimaliz´al´asa v´egett oldott´ak meg a r´ezsz´al vagy r´ezsz´alak k¨orbetekercsel´es´evel [5]. Hiszen ilyenkor a h´ ur merevs´ege kisebb, mintha egy ugyanolyan t¨omeg˝ u, de t¨om¨or h´ urt haszn´aln´anak a gy´art´ok. A nagyobb s´ ulyra az´ert van sz¨ uks´eg, mert a h´ urok k¨ozti m´eretk¨ ul¨onbs´egek korl´atozottak, valamint a fesz´ıt´es m´ert´ek´eben sem lehetnek nagy k¨ ul¨onbs´egek, hogy a t˝ok´et ne ´erje kiegyens´ ulyozatlan er˝ohat´as [5]. A m´eretk¨ ul¨onbs´egek persze egy koncertzongor´an´al j´oval nagyobbak lehetnek egy pian´ın´ohoz k´epest, ´es ez meg is nyilv´anul az inharmonicit´as m´ert´ek´eben [12]. A hetvenes ´evek ´ota tudjuk, hogy a zongorahang spektrum´aban l´eteznek egy´eb felhangok is, melyek inharmonicit´asa k¨or¨ ulbel¨ ul a negyede az el˝obb ismertetetteknek [27]. Ezeket nevezi ma m´ar a szakirodalom ,,fantom felhangoknak” (phantom partials) [27]. [3] szerint ezeket a h´ urban keletkez˝o longitudin´alis rezg´esek seg´ıts´eg´evel lehet magyar´azni, s˝ot a m´er´esekkel j´ol egyez˝o modellt is ismertet. [2]-ben azt vizsg´alj´ak a szerz˝ok, hogy az ilyen rezg´esek hat´asa mennyire hallhat´o, ´es szint´en arra jutottak, hogy a m´elyebb hangok eset´en ´erv´enyes¨ ulnek ilyen ´ertelemben a longitudin´alis vibr´aci´ok. Term´eszetesen figyelembe lehet venni a csillap´ıt´ast is, mely h´ ur eset´en okozhat a k¨orbefoglal´o leveg˝o viszkozit´asa, a h´ ur bels˝o fesz¨ ults´ege, ´es a hidakon leadott energia ´ elveszt´ese [12]. Erdekes, hogy a zongorahang nem folyamatos u ¨temben hal el, hanem k´et r´eszre oszthat´o: ,,[...] a korai szakaszban gyorsabban hal el, mint a k´es˝oiben”[1]. Szok´as ezt ,,dupla elhal´asnak” (two-stage decay) nevezni[1]. Csak a legegyszer˝ ubb (3.1) fizikai modellt vizsg´aljuk a tov´abbiakban. 10
3.2.
A rezon´ ans modellj´ enek fel´ır´ asa
3.2.1.
A rezon´ ans feladata ´ es fel´ ep´ıt´ ese
A rezon´ans a h´ urok ´es a h´ıd mozg´asi energi´aj´anak egy r´esz´et akusztikus energi´av´a alak´ıtja [12]. Azaz a h´ urok hangj´at feler˝os´ıti. Szinte minden h´ uros hangszer eset´en ez a f˝o feladata. Fontos megeml´ıteni, hogy sok h´ uros hangszer eset´en ezt nem csak egy lap v´egzi, hanem egy ,,lyukas doboz”. Ilyenek p´eld´aul a (klasszikus ´ertelemben vett) von´osok, a git´ar, a h´arfa. Az ilyen t´ıpus´ u hangszertestek modellez´ese t´ ulmutat jelen munka keretein. A zongora eset´en is el˝ofordul´o esettel foglalkozom, amikor csak egy rezon´al´o lap van. Ez a modell a git´ar k¨ozel´ıt´es´ere is haszn´alhat´o ha megel´egsz¨ unk a magasabb frekvenci´ak kelt´es´evel. Ezt a meg´allap´ıt´asomat arra a t´enyre alapozom, hogy eme hangszer eset´en az eml´ıtett frekvenci´akat l´enyeg´eben az el˝olap sug´arozza [12]. A bonyolultabb eset ir´ant ´erdekl˝od˝o olvas´onak javaslom a k¨ovetkez˝o forr´asokat: [12, 9, 10]. A h´ urok ´altal kifejtett er˝o a h´ıdon kereszt¨ ul ´eri a rezon´anst. Ez az o¨sszekapcsol´as tov´abb´ıtja a rezg´eseket de azt is okozza, hogy a h´ urokban lev˝o fesz¨ ults´egb˝ol ad´od´o ´alland´o er˝o is hat a rezon´ansra. Ez a h´ urok a´ltali a´lland´o terhel´es a zongor´ak, pian´ın´ok eset´eben olyan er˝os, hogy a gyeng´ebb min˝os´eg˝ uekn´el a rezon´ans ak´ar el is deform´al´odhat [12]. A zongora ill. pian´ın´o rezon´ans´anak fel´ep´ıt´es´et [12] alapj´an ismertetem. A zongora rezon´ansa nem teljesen lapos, a h´ıd alatt van egy 1–2 mm-es korona, valamint a sz´eleken ´altal´aban elv´ekonyodik. Anyag´at tekintve p´arhuzamosan o¨sszeragasztott lucfeny˝o cs´ıkokb´ol van. A cs´ıkokra mer˝oleges ir´anyban a rezon´ans h´atoldal´ara bord´ak vannak er˝os´ıtve. Ezek a farostok kisebb keresztir´any´ u merevs´eg´et hivatottak p´otolni. Egy viszonylag friss kutat´as m´er´esei r´amutatnak arra, hogy ezt a terhel´es ´ is seg´ıti [24]. Erdekes k´erd´es, hogy hogyan van r¨ogz´ıtve a rezon´ans. Ez [12]-ben nincs r´eszletezve. A rezon´ans h´ urok fel˝oli r´eszein szeg´ely van, melyek seg´ıts´eg´evel a t´ uloldalon lev˝o vastag fa zongoratesthez csavarozz´ak. A csavarok egym´as k¨oz¨otti t´avols´aga az ´en pian´ın´om eset´en kb. 20 cm. [22] eml´ıti m´eg, hogy a rezon´anst el˝o szokt´ak fesz´ıteni oly m´odon, hogy a bord´ak ragaszt´asi fel¨ ulete kicsit hajl´ıtott.
3.2.2.
Eddigi eredm´ enyek a fizikai modellek kutat´ as´ aban
A rezon´ans modellj´et a h´ ur ´es a kalap´acs modellj´ehez viszony´ıtva keveset vizsg´alt´ak (hasonl´o meg´allap´ıt´ast tal´alunk [14, 15]-ben is, a helyzet az´ota l´enyeg´eben nem v´altozott), b´ar az ut´obbi t´ız ´evben ennek a kutat´asa is felgyorsult. Az 1998-as kiad´as´ u [12] az alkatr´esz mod´alis anal´ızis´evel kapcsolatban l´enyeg´eben csak n´egy
11
cikkre hivatkozik, melyek egyike sem tud felmutatni m´er´esekkel j´ol egyez˝o numerikus modellt. Ismereteim szerint az´ota sem siker¨ ult ilyet el˝oa´ll´ıtani annak ellen´ere, hogy az ut´obbi id˝oszak kutat´asait eln´ezve egy´ertelm˝ u a fejl˝od´es e t´eren is. Kezdetben a rezon´anst o¨nmag´aban (de a h´ıddal egy¨ utt) vizsg´alt´ak [11, 29, 13, 6, 14, 15], k´es˝obb elkezdt´ek vizsg´alni a h´ urok terhel´es´evel is sz´amolva [23, 24]. A legut´obbi kutat´asok arra is keresik a v´alaszt, hogy milyen hat´assal van a saj´atrezg´esek alakul´as´ara a rezon´ans el˝ofesz´ıt´ese [22]. Az o¨nmag´aban val´o vizsg´alatkor legt¨obbsz¨or az eredeti hely´er˝ol kiemelt´ek az alkatr´eszt, m´ıg a terhel´eseket figyelembe v´eve a hangszert szinte ´erintetlen¨ ul vizsg´alt´ak, az el˝ofesz´ıt´est viszont egy egyszer˝ us´ıtett hangszeren [22].
3.2.3.
Az egydimenzi´ os modell
Fizikailag a lemez egydimenzi´os megfelel˝oje a gerenda (bar) vagy r´ ud (rod). Felfog´ hat´o u ´gy mint egy merevs´eggel rendelkez˝o h´ ur. ,,Altal´ anoss´agban a fesz´ıt´es, mint visszah´ uz´o er˝o fontosabb a merevs´egn´el a h´ ur eset´en ´es a merevs´eg fontosabb a gerenda eset´en [...].” [26]. Mi most a gerend´ara ´es a lemezre is u ´gy tekint¨ unk, amiben nincs fesz´ıt´es, azaz az egyetlen er˝o, ami az egyens´ ulyi a´llapotba viszi vissza, a merevs´egb˝ol ad´odik. Ha egy´eb jelens´egeket (mint pl. a longitudin´alis mozg´asok, csillap´ıt´as) nem vesz¨ unk figyelembe, akkor a mi eset¨ unkben a gerenda mozg´asegyenlete a k¨ovetkez˝o alakot o¨lti [26, 12]: ρ
4 ∂ 2u 2∂ u + Eκ =f dt2 dx4
(3.3)
ahol ρ a s˝ ur˝ us´eget E a Young-f´ele modulust κ a forg´as sugar´at jel¨oli, u : Ω × [0, T ] → R a keresett megold´as az Ω × [0, T ] tartom´anyon, T ∈ R, Ω tartom´any, f : Ω×[0, T ] → R adott. Az Ω tartom´any egy z´art [a, b] intervallum (a, b ∈ R, a < b). Az egyenlet levezet´ese megtal´alhat´o a [26, 12] m˝ uvekben. Nem homog´en esetben a jobboldal itt is egy k¨ uls˝o er˝ot tud bevinni az egyenletbe. A peremfelt´etelek azonosak a (3.6)-ban le´ırtakkal az ottani megfontol´asok alapj´an. Term´eszetesen az itteni Ω tartom´anyt kell figyelembe venni. A kezdeti´ert´ekek a h´ uregyenlethez hasonl´oan a k¨ovetkez˝ok: u(x, 0),
∂u (x, 0) (x ∈ [a, b]) ∂x
(3.4)
Azaz kezdetben adottak a kit´er´esek ´es a sebess´egek. Felmer¨ ul a k´erd´es, hogy ´erdemes-e ezzel a modellel foglalkozni. Ez a modell egyr´eszt seg´ıtheti a k´etdimenzi´os megold´asra jut´ast, m´asr´eszt ki tudja fejezni a h´ urok egyfajta ¨osszekapcsol´as´at, s˝ot lehet a h´ıd modellje is. Ha a fel´ırt gerend´ara 12
r´ak¨ot¨ unk h´ urokat, akkor az olyan, mintha nem venn´enk figyelembe a rezon´anst, de a hidat igen. Ez egyr´eszt hasznos lehet olyan esetben, amikor csak a h´ ur mozg´asa ´erdekel minket, de valamennyire re´alis felf¨ uggeszt´essel. M´asr´eszt seg´ıts´eget adhat olyan k´erd´esek megv´alaszol´as´ara, amelyek a h´ ur h´ıdra val´o hat´as´at vizsg´alj´ak. V´eg¨ ul n´ezetem szerint viszonylag j´ol le lehet ´ırni vele az unisono h´ urok o¨sszekapcsol´as´at. Ezek a h´ urok igen k¨ozel vannak egym´ashoz (n´eh´any millim´eter). A fenti esetek nem ´erdektelenek, hisz a szakirodalom alapj´an kider¨ ul, hogy a kutat´o t´arsadalom folyamatosan vizsg´alja ezeket.
3.2.4.
A k´ etdimenzi´ os modell
A legegyszer˝ ubb k´etdimenzi´os modell a rezon´ansra a v´ekony izotr´op lemez. [24] alapj´an ez a modell kieg´esz´ıtve a hidakkal ´es terhel´essel eg´esz j´o egyez´est mutat legal´abb az els˝o h´arom saj´atfrekvenci´at tekintve, a saj´atf¨ uggv´enyek viszont m´ar kev´esb´e egyeznek. A v´ekony izotr´op lemez mozg´asegyenlete [12, 26]: ρs
∂ 2u Eq 3 + ∆2 u = f dt2 12(1 − ν 2 )
(3.5)
ahol ρ a s˝ ur˝ us´eg, E a Young-f´ele modulus, s a lemez vastags´aga, ν a Poisson ar´any, ∆2 a biharmonikus oper´ator (∆ a Laplace oper´ator, ∆ = ∇2 ), u : Ω × [0, T ] → R a keresett megold´as Ω×[0, T ] tartom´anyon, T ∈ R, Ω ⊂ R2 tartom´any, f : Ω×[0, T ] → R adott. Az egyenlet levezet´ese megtal´alhat´o [26]-ban. A biharmonikus oper´ator ki´ırva: ∆2 u =
∂ 4u ∂ 4u ∂ 4u + 2 + ∂x4 ∂x2 ∂y 2 ∂y 4
Megjegyzem tov´abb´a, hogy a szakirodalom s-t h-val szokta jel¨olni, de mivel ´en azt a diszkretiz´aci´o l´ep´est´avols´ag´ara haszn´alom, ez´ert m´as karaktert v´alasztottam. Arra val´oban lehet rem´eny¨ unk, hogy a nagy hull´amhossz´ u rezg´eseket viszonylag j´ol visszaadja a modell, hiszen ezek j´oval nagyobbak a bord´ak t´avols´ag´an´al, a kisebb hull´amhosszok eset´en viszont elker¨ ulhetetlen az ortotr´op eset ´es m´eg a bord´ak figyelembe v´etele is, azaz inhomog´en lesz a lemez [13]. Mindez abb´ol ad´odik, hogy az igazi rezon´ans a nagy hull´amhossz´ u rezg´esek eset´en merevebb lemezk´ent viselkedik, mint a bord´ak egym´as k¨ozti t´avols´ag´an´al kisebb hull´amhossz´ uak eset´en [13]. Nem volt sz´o m´eg a peremfelt´etelekr˝ol. Ismereteim szerint ezt a k´erd´est a tudom´any m´eg nem v´alaszolta meg. Egy a [12]-ben is ismertetett kutat´as m´er´esei a pian´ın´o eset´eben egyfajta a´tmeneti esetet sugallnak a befogott ´es az al´at´amasztott sz´eleknek megfelel˝o peremfelt´etelek k¨oz¨ott. A legt¨obb tanulm´any azonban befogott 13
sz´eleket felt´etelez, mint p´eld´aul [11, 13, 24]. A saj´at pian´ın´omat eln´ezve v´elem´enyem szerint ez a feltev´es nem ´all t´ ul messze a val´os´agt´ol. El´eg ha csak a szeg´elyekre gondolunk, igaz vannak r´eszek, ahol nincs szeg´ely. A befogott sz´eleknek megfelel˝o peremfelt´etelek a k¨ovetkez˝ok: u|Γ = 0 ∂u |Γ = 0 dt
(3.6)
ahol Γ az Ω tartom´any perem´et jel¨oli. Most csak a n´egyzetes tartom´anyt vizsg´alom: Ω = [a, b] × [a, b] (a, b ∈ R). Fontos vizsg´alni a rezon´ans saj´atfrekvenci´ait ´es saj´atf¨ uggv´enyeit. A k´erd´es a ∆2 oper´ator saj´at´ert´ek-probl´em´aja. Val´oban, ha (3.5) megold´as´at u(x, y, t) = Z(x, y)eiωt alakban keress¨ uk f ≡ 0 eset´en, akkor az −
12ρ(1 − ν 2 )ω 2 Z(x, y) + ∆2 Z(x, y) = 0 2 Es
egyenl˝os´egre jutunk [12]. Ezt a´trendezve kapjuk a ∆2 Z(x, y) =
12ρ(1 − ν 2 )ω 2 Z(x, y) Es2
(3.7)
azonoss´agot. A jobboldalon Z egy¨ utthat´oja konstans. Sajnos az ´altalunk v´alasztott tartom´any ´es peremfelt´etelek eset´en nem l´etezik a feljebbi saj´at´ert´ekprobl´ema z´art alakban fel´ırhat´o megold´asa. Viszont ismertek k¨ozel´ıt´esek, becsl´esek mind az ´ert´ekekre, mind a f¨ uggv´enyekre. Ilyeneket k¨oz¨ol [12, 7]. Bevezetem m´eg a szakirodalom ´altal gyakran haszn´alt D :=
Es2 12ρ(1 − ν 2 )
jel¨ol´est.
14
4. fejezet Az ¨ osszetev˝ ok ¨ osszekapcsol´ as´ anak modelljei 4.1.
A h´ ur ´ es a rezon´ ans k¨ ozvetlen ¨ osszekapcsol´ asa
Ebben a pontban egy olyan modellt ismertetek, mely a h´ ur v´egpontj´at a rezon´ans egy pontj´aval kapcsolja o¨ssze az´altal, hogy a k´et pont gyorsul´as´at egyenl˝ov´e teszi.
4.1.1.
Az egydimenzi´ os eset
Tekints¨ uk az (3.1) ´es (3.3) mozg´asegyenleteket. Mindkett˝oben a jobboldalt tekints¨ uk azonosan null´anak. Az azonos nev˝ u v´altoz´okat, feloszt´asokat indexelj¨ uk a h´ ur eset´en st-vel, a gerenda eset´en r-rel. A gerenda egy adott xr,i pontj´ahoz fogjuk kapcsolni a h´ ur xst,L v´egpontj´at. Tegy¨ uk fel, hogy ∂ 2 ur ∂ 2 ust (x , t) = (xst,L , t) (t ∈ [0, T ]). r,i ∂t2 ∂t2
(4.1)
Ebb˝ol az o¨sszef¨ ugg´esb˝ol ´es a mozg´asegyenletekb˝ol kapjuk, hogy Eκ2 ∂ 4 ur Θ ∂ 2 ust (x , t) = − (xst,L , t) (t ∈ [0, T ]). r,i ρ ∂x4 µ ∂x2
4.1.2.
(4.2)
A k´ etdimenzi´ os eset
A h´ ur ´es a lemez ¨osszekapcsol´asa ´erdek´eben tekints¨ uk az (3.1) ´es (3.5) mozg´asegyenleteket. Mindkett˝oben a jobboldalt tekints¨ uk azonosan null´anak. Az egydimenzi´os esethez hasonl´oan fogunk elj´arni. Az azonos nev˝ u v´altoz´okat, feloszt´asokat indexelj¨ uk a h´ ur eset´en tov´abbra is st-vel, a lemez eset´en pl-lel. A lemez egy adott (xpl,i , ypl,j ) pontj´ahoz fogjuk kapcsolni a h´ ur xst,L v´egpontj´at. Tegy¨ uk fel, hogy ∂ 2 upl ∂ 2 us t (x , t) = (xst,L , t) (t ∈ [0, T ]). pl,i ∂t2 ∂t2 15
(4.3)
Ebb˝ol az o¨sszef¨ ugg´esb˝ol ´es a mozg´asegyenletekb˝ol kapjuk az Θ ∂ 2 ust (xst,L , t) (t ∈ [0, T ]). µ ∂x2
−D∆2 upl (xpl,i , ypl,j , t) =
(4.4)
azonoss´agot.
4.2.
A h´ ur ´ es a rezon´ ans o ¨sszekapcsol´ asa egy egyszer˝ u h´ıdon kereszt¨ ul
Ebben a r´eszben egy a [4]-ben k¨oz¨olt m´odszert alkalmazom a mi modell¨ unkre. [4] a git´ar h´ıdj´at u ´gy modellezi, hogy az el¨ uls˝o lap a h´ıd a´ltal fedett r´eszen integr´alt sz´am´ıt. Ez az integr´al azt fejezi ki, hogy a h´ıd egyfajta ´atlagol´ast v´egez. A modellben a h´ıd t¨omege nem jelenik meg. Csak az egydimenzi´os esetet dolgozom ki. Tekints¨ uk az (3.3), (3.1) mozg´asegyenleteket, ´es haszn´aljuk itt is a 4.1.1-beli jel¨ol´eseket. Vezess¨ uk be a G s´ ulyf¨ uggv´enyt, mely a gerenda azon r´esz´et jel¨oli ki, amely f¨ol¨ott szeretn´enk integr´alni [4]-hez hasonl´oan. b
Z G : [a, b] → R,
G(x)dx = 1 a
Mi a legegyszer˝ ubb G(x) =
0
ha x ∈ [a, v) 1
ha x ∈ [v, w]
w−v 0
(4.5)
ha x ∈ (w, b]
esettel foglalkozunk, ahol v, w ∈ [a, b] a h´ıd (´es ´ıgy az ,,´atlagol´as”) k´et sz´ele (v < w). A h´ ur a gerend´ara a gerenda mozg´asegyenlet´enek jobboldal´an kereszt¨ ul hat, ´ıgy az a k¨ovetkez˝ok´eppen n´ez ki [4]: fr (x, t) = −Θ
∂ust (L, t)G(x). ∂xst
(4.6)
G elosztja a h´ ur er˝ohat´as´at a rezon´ator h´ıddal fedett r´esz´en [4]. A gerenda a h´ urra a kit´er´es´en kereszt¨ ul hat, de a h´ıd k¨ozvet´ıt´es´evel [4]: Z ust (L, t) =
b
G(xr )ur (xr , t)dxr . a
[4] bizony´ıtja az ilyen ¨osszekapcsol´as energiamegmarad´as´at.
16
(4.7)
4.3.
Az egyszer˝ u h´ıdon kereszt¨ ul val´ o ¨ osszekapcsol´ as egy v´ altozata
Haszn´aljuk az el˝oz˝o pont jel¨ol´eseit. Az el˝oz˝o pontok alapj´an ad´odik az ¨otlet, hogy a gerenda egy r´esze sebess´egei f¨ol¨otti normaliz´alt integr´alt tegy¨ uk egyenl˝ov´e a h´ ur v´egpontj´anak sebess´eg´evel. Azaz bizonyos ´ertelemben o¨tv¨ozz¨ uk a 4.1.1-ben ´es a 4.2ben bemutatott m´odszereket. Form´alisan megfogalmazva: Z b ∂ 2 ur (xr , t) ∂ 2 ust G(xr ) dxr = (t, L). ∂t2 ∂t2 a
(4.8)
Mind k´et oldalon a m´asodrend˝ u deriv´altakat kifejezve a nekik megfelel˝o mozg´asegyenletekb˝ol kapjuk a k¨ovetkez˝o egyenl˝os´eget: Z
G(xr ) a
4.4.
2
4
b
fr − Eκ2 ∂∂xu4r r
ρ
dxr =
fst + Θ ∂xu2 st st
µ
.
(4.9)
A gerjeszt´ es modellez´ ese
A szakirodalom h´arom gerjeszt´est szokott megk¨ ul¨onb¨oztetni a h´ ur eset´eben: penget´es, (kalap´accsal val´o) u ¨t´es, von´as. Pengetni p´eld´aul git´art, h´arf´at, banjot szokt´ak, de a csembal´o eset´eben is err˝ol besz´elhet¨ unk, csak ott nem k¨ozvetlen¨ ul a j´at´ekos penget, hanem egy k¨ozvet´ıt˝o mechanika seg´ıts´eg´evel egy penget˝o (eredetileg l´ udtoll volt). A zongora eset´eben egy kalap´acs u ¨ti meg a h´ urt. A kalap´acsot a j´at´ekos egy bonyolult mechanik´an kereszt¨ ul tudja leng´esbe hozni. Von´asr´ol a von´osok eset´eben besz´el¨ unk, b´ar a kort´ars zen´eben sok m´as eszk¨ozt is szoktak vonni (p´eld´aul f˝ ur´eszt, cint´any´ert, git´ar h´ urj´at). A von´as modellez´es´evel annak bonyolults´aga miatt nem foglalkozom. Ha az olvas´o ezt az esetet meg szeretn´e ismerni, javaslom a [12] o¨sszefoglal´o m˝ uvet elindul´ask´ent.
4.4.1.
Kezdeti ´ ert´ ekkel val´ o gerjeszt´ es
A legegyszer˝ ubb m´odja a h´ ur gerjeszt´es´enek, hogy a kezdeti ´ert´ekeket a´ll´ıtjuk be u ´gy, hogy a nek¨ unk megfelel˝o esetet k¨ozel´ıts¨ uk. [12] alapj´an k¨ozl¨om a vizsg´alatunk t´argy´aul szolg´al´o k´et esetet: Penget´est akkor k¨ozel´ıt¨ unk, ha a kezdeti kit´er´est megv´altoztatjuk az egyens´ ulyi a´llapothoz k´epest, a kezdeti sebess´egek null´ak. A kezdeti kit´er´est u ´gy kell meghat´arozni, hogy a h´ ur kezdeti ´es egyens´ ulyi kit´er´ese egy h´aromsz¨oget alkosson. A penget´es hely´et˝ol is f¨ ugg hogy mely felhangok fognak sz´olni. Egy ilyen esetet mutat be a 4.1 ´abra. 17
4.1. ´abra. Gerend´ara kapcsolt k´et h´ ur k¨oz¨ ul az egyik penget´essel gerjesztve. A piros pontok a kezdeti kit´er´est mutatj´ak. Ahol nincsenek, ott a kezdeti kit´er´es azonosan nulla volt. Diszkretiz´aci´os param´eterek a k¨ovetkez˝o fejezetben bevezetett jel¨ol´esekkel: n = 20, m = 800. Az ´abra a k = 5080 id˝or´eteget mutatja. A pengetett h´ ur konstansai Θ = 1, µ = 1, a m´asik h´ ur´e Θ = 1.58, µ = 1, a gerend´a´e Eκ2 = 0.5, ρ = 1. T = 1. A kapcsol´od´asi pontok rendre x6 , x12 . Kalap´accsal u ¨t´est akkor k¨ozel´ıt¨ unk ha a kezdeti kit´er´es az egyens´ ulyi a´llapot, viszont a sebess´eg egy pontban nem nulla [12, 1].
4.4.2.
A jobboldalon kereszt¨ ul t¨ ort´ en˝ o gerjeszt´ es
Jobb k¨ozel´ıt´est lehet el´erni, ha (3.1)-nek az f jobboldal´an egy, a gerjeszt˝o eszk¨oz (p´eld´aul kalap´acs, j´at´ekos ujja) a´ltal a h´ urra kifejtett er˝ot alkalmazunk. A git´ar ujjal val´o penget´es´enek k¨ ul¨onb¨oz˝o modelljeit megtal´alhatjuk [8, 4] m˝ uvekben. A kalap´accsal val´o u ¨t´es t¨obb modellj´et megtal´alhatjuk [1, 27, 16]-ban. A penget´es eset´en egyszer˝ ubb ezt az er˝ot le´ır´o f¨ uggv´enyeket tal´alunk, az u ¨t´es eset´en implicit formul´akat tal´alunk, melyek sz´am´ıt´asa numerikus probl´em´akat (instabilit´asit is) felvet. Az eml´ıtett numerikus probl´em´akr´ol [1] ´ır.
18
5. fejezet Numerikus megval´ os´ıt´ asok A megval´os´ıt´asok sor´an a legegyszer˝ ubb, explicit v´egesdifferencia s´em´akkal dolgozom. H´atr´anyuk, hogy a stabilit´asuk sok id˝or´eteget k´ıv´an, el˝ony¨ uk viszont az, hogy p´arhuzamos´ıt´asuk term´eszetes m´odon ad´odik, b´ar a megval´os´ıt´as sor´an ezt nem haszn´altam ki. Tov´abbi, m´eg fontosabb el˝ony¨ uk, hogy a nagy frekvenci´aj´ u megold´asokat jobban visszaadj´ak, mint az implicit m´odszerek. Az eg´esz fejezetre igaz, hogy a tesztek sor´an igyekeztem min´el egyszer˝ ubb f¨ uggv´enyeket haszn´alni, ez´ert ahol nem jelzem k¨ ul¨on, ott a mozg´asegyenletek konstansait u ´gy hat´aroztam meg, hogy a parci´alis deriv´altak egy¨ utthat´oi az 1 ´ert´eket kapj´ak. A tesztek eredm´enyeit bemutat´o t´abl´azatokban a ,,Hiba”, maxim´alis hib´at jelent, a ,,Sebess´eg” pedig konvergenciasebess´eget, mely az el˝oz˝o, ´es az aktu´alis sor hib´aj´anak ar´anya. Az´ert ez a konvergenciasebess´eg, mert az egy sorhoz tartoz´o approxim´aci´o egy dimenzi´o ir´any´aba val´o t´erbeli feloszt´asa oszt´opontjainak sz´ama n´alunk a fele az alatta lev˝o sor eset´ehez k´epest.
5.1. 5.1.1.
A h´ ur egy egyszer˝ u modellj´ enek megval´ os´ıt´ asa A v´ egesdifferencia m´ odszer fel´ır´ asa
A h´ ur eset´eben tekints¨ uk a (3.1) mozg´asegyenletet a hozz´atartoz´o jel¨ol´esekkel. Az [a, b] ill. [0, T ] intervallumot osszuk fel n ∈ N ill. m ∈ N r´eszre ekvidiszt´ans m´odon. Jel¨olj¨ uk x0 , x1 , ..., xn -nel ill. t0 , t1 , ..., tm -nel az oszt´opontokat ´es h-val ill. τ -val a szomsz´edos pontok k¨ozti t´avols´agot. Vil´agos, hogy h = (b − a)/n ill. τ = T /m. A h´ ur eset´eben csak az a = 0, b = L speci´alis esettel foglalkozom. Ekkor h = L/n. uki -val jel¨ol¨om u(xi , tk )-t illetve az azt approxim´al´o v´altoz´ot, att´ol f¨ ugg˝oen, hogy hol haszn´alom. A mozg´asegyenlet alapj´an a v´egesdifferencia s´em´ahoz sz¨ uks´eg van a m´asodrend˝ u
19
deriv´alt szok´asos approxim´aci´oj´ara: ukt¯t,i,τ :=
1 k−1 (u − 2uki + uk+1 ) i τ2 i
(5.1)
A l´ep´est´avols´ag (itt most τ ) jel¨ol´es´et elhagyom, kiv´eve, ha nem egy´ertelm˝ u. [17] alapj´an ez a differencias´ema m´asodrend˝ u approxim´aci´ot ad, felt´eve, hogy u(xi , t) ∈ C 4 [0, T ]. A fentieket felhaszn´alva kapjuk a (3.1) diszkretiz´aci´oj´at: ukt¯t,i −
Θ k ux¯x,i = fij µ
(i = 1, 2, . . . , n − 1, k = 1, 2, . . . , n − 1) Ebb˝ol kifejezve uk+1 -et az i uk+1 = 2uki − uk−1 + i i
Θτ 2 k τ2 ux¯x,i + fik µ µ
(i = 1, 2, . . . , n − 1, k = 1, 2, . . . , n − 1) m´asodrend˝ u explicit s´em´ara jutunk. Hi´anyzik m´eg a k = 0 eset, azaz a k = 1 id˝or´eteg kisz´am´ıt´as´anak s´em´aja. Ez [18]-ban r´eszletesen le van vezetve, ´en csak v´azlatosan k¨ozl¨om az ottani levezet´est. El˝osz¨or is azt kell l´atni, hogy ∂u 1 u1i − u0i = (xi , ) + O(τ 2 ) (i = 0, 1, . . . , n). τ ∂t 2 Ez ad´odik az u-nak az (xi , 1/2) k¨or¨ uli t ir´any´ u sorfejt´es´eb˝ol [17]-ben r´eszletezett m´odon. A m´asik o¨tlet ugyancsak Taylor-sor alapj´an ad´odik: ∂u 1 ∂u τ ∂ 2u (xi , ) = (xi , 0) + (xi , 0) + O(τ 2 ) (i = 0, 1, . . . , n). 2 ∂t 2 ∂t 2 ∂t A k´et azonoss´agb´ol ad´odik, hogy az el˝oz˝o baloldala megegyezik az ut´obbi jobboldal´aval. A m´asodrend˝ u id˝oderiv´alt viszont kifejezhet˝o a mozg´asegyenletb˝ol, ´es ´ıgy approxim´alhat´o is, viszont a hibatag O(τ h+τ 2 ) lesz. V´eg¨ ul csak arra kell hivatkozni, hogy τ h ≤ (τ 2 + h2 )/2 ´es ezzel bel´attuk a k¨ovetkez˝o s´ema m´asodrend˝ us´eg´et [18]: u1i = u0i + τ
∂u τ2 (xi , 0) + (Θu0x¯x,i + fi0 ) (i = 1, 2, . . . , n − 1). ∂t 2µ
3.1 peremfelt´etelei, kezdeti´ert´ekei, jobboldala form´alisan most azt jelentik, hogy adottak az u0i (i = 0, 1, . . . , n) ∂u u01i := (i = 0, 1, . . . , n) ∂t uk0 , ukn (k = 1, 2, . . . , m) ´ert´ekek. 20
n
m
Hiba
10
200
0.0135868
20
800
0.0034443
3.9447202
40
3200
0.000864103
3.98598312
80
12800
0.000216213
3.996535823
51200
−5
3.999130676
160
Sebess´eg
5.4065·10
5.1. t´abl´azat. A h´ ur megval´os´ıt´as´anak hib´aja (µ = 1, Θ = 1). n
m
Hiba
Sebess´eg
10
200
0.0195229
20
800
0.00494192
3.950468644
40
3200
0.00123931
3.987638283
80
12800
0.000310067
3.996910345
160
51200
7.75317·10−5
3.999228703
5.2. t´abl´azat. A h´ ur megval´os´ıt´as´anak hib´aja (µ = 1, Θ = 1.5).
5.1.2.
Stabilit´ as
A m´odszer stabilit´as´anak vizsg´alata a Neumann-f´ele m´odszerrel [18]-ban alaposan r´eszletezve van. Az ottani eredm´eny szerint a v´eges intervallumon vett stabilit´as el´egs´eges felt´etele a k¨ovetkez˝o: |c|τ <1 h ahol c =
5.1.3.
q
(5.2)
Θ . µ
Tesztek
T¨obb tesztf¨ uggv´enyre kipr´ob´altam a m´odszert, ezek k¨oz¨ ul egyet mutatok be, de a feladatot k´etf´ele m´odon param´eterezem. Az els˝o param´eterez´es szerint homog´en a m´asodik szerint inhomog´en egyenletre jutunk. A tesztf¨ uggv´eny: u = cos t cos x (x, t ∈ [0..2π]). A megval´os´ıt´asok hib´aj´at az 5.1, 5.2 t´abl´azatok tartalmazz´ak. J´ol l´athat´o a m´asodrend˝ u konvergencia.
21
5.2.
Az egydimenzi´ os rezon´ ans egyszer˝ u modellj´ enek megval´ os´ıt´ asa
5.2.1.
A v´ egesdifferencia m´ odszer fel´ır´ asa
Tekints¨ uk [a, b] 5.1-ben tal´alhat´o feloszt´as´at ´es annak jel¨ol´eseit. A gerenda (3.3) mozg´asegyenlet´et ezen a feloszt´ason diszkretiz´alva sz¨ uks´eg¨ unk van a m´asodrend˝ u differenci´alh´anyados (5.1) approxim´aci´oj´an t´ ul a negyedrend˝ u differenci´alh´anyados approxim´aci´oj´ara is. Ezt a [17]-ben k¨oz¨olt Taylor-sorfejt´esen alapul´o m´odszerrel nem neh´ez levezetni, a hivatkozott m˝ u ´altal k¨oz¨olt eredm´eny: ukx¯x¯xx,i,h :=
1 k (u − 4uki−1 + 6uki − 4uki+1 + uki+2 ). h4 i−2
(5.3)
Enn´el a jel¨ol´esn´el is igaz, hogy ha a l´ep´est´avols´agot nem sz¨ uks´eges ki´ırni, akkor elhagyom. Az eml´ıtett levezet´esb˝ol ad´odik, hogy az approxim´aci´o m´asodrend˝ uu∈ C 6 ([0, L]) eset´en [17]. A gerenda mozg´asegyenlet´enek diszkretiz´aci´oja ezek ut´an a k¨ovetkez˝ok´eppen ´ırhat´o fel: ρukt¯t,i + Eκ2 ukx¯x¯xx,i = fik (i = 2, . . . , n − 2,
(5.4)
k = 1, . . . , m − 1).
Ha kifejezz¨ uk uk+1 -et, akkor a k¨ovetkez˝o explicit s´em´at kapjuk: i uk+1 = 2uki − uk−1 − i i
Eκ2 τ 2 k τ 2 fik ux¯x¯xx,i + ρ ρ
(i = 2, . . . , n − 2,
k = 1, . . . , m − 1).
(5.5)
A k = 1-hez tartoz´o id˝or´eteg kisz´amol´as´ahoz a k¨ovetkez˝o s´em´at alkalmazzuk: u1i
=
u0i
∂u Eκ2 τ 2 0 τ 2 fi0 + τ (xi , 0) − u + ∂t 2ρ x¯x¯xx,i 2ρ
(5.6)
(i = 2, . . . , n − 2). Ezt a s´em´at a hull´amegyenlet´ehez hasonl´oan lehet levezetni, a m´asodrend˝ us´ege is hasonl´oan ad´odik. Ez ut´obbi esetben csak azt a t´enyt kell alkalmazni, hogy a negyedrend˝ u deriv´altat approxim´al´o s´ema (megfelel˝o differenci´alhat´os´ag eset´en) ugyan´ ugy m´asodrend˝ u, mint a hull´amegyenletn´el alkalmazott m´asodrend˝ u helyderiv´altat approxim´al´o. A kezdeti´ert´ekekk´ent ´es a peremfelt´etelekk´ent (3.4) ´es (3.6) alapj´an a k¨ovetkez˝o ´ert´ekek adottak: u0i , uk0 , ukn ,
∂u (xi , 0) (i = 0, 1, . . . , n) ∂t
∂u ∂u (x0 , tk ) = 0, (xn , tk ) = 0 (k = 0, 1 . . . , m). ∂x ∂x 22
(5.7)
Adott m´eg fik (i = 0, 1, . . . , n, k = 0, 1 . . . , m). A kapott k´epletekben a fut´o indexeket megvizsg´alva kider¨ ul, hogy az i = 1, n − 1 helyek kisz´am´ıt´as´at nem ´ırj´ak le. Ez nem v´eletlen, hisz a negyedrend˝ u deriv´alt o¨tpontos s´em´aj´at ezekre a pontokra nem tudjuk alkalmazni. ´Igy a peremen bel¨ uli els˝o helyeket az adott id˝or´etegben m´ar kisz´amolt k¨ornyez˝o pontok seg´ıts´eg´evel kell k¨ozel´ıteni. Vegy¨ uk p´eld´aul az u1 ´ert´eket. u1 ´es u2 -nek az x0 k¨or¨ uli Taylor-sorba fejt´es´et felhaszn´alva kapjuk a k¨ovetkez˝o azonoss´agot: 1 3 1 ∂u 1 2 ∂ 3u (− u0 + 2u1 − u2 ) = (0) − h (0) + O(h3 ) 3 h 2 2 ∂x 3 ∂x Mivel
∂u (0) ∂x
(5.8)
a mi eset¨ unkben 0, ez´ert u1 -re rendezve az al´abbi approxim´aci´ot kapjuk: 1 u1 = (3u0 + u2 ) + O(h3 ) 4
(5.9)
(A gerenda m´asik v´eg´ere az approxim´aci´o az el˝obbi k´epletre szimmetrikus.) Teh´at a k. id˝or´etegben az al´abbi k¨ozel´ıt´est v´egezz¨ uk: 1 uk1 := (3uk0 + uk2 ) 4
(5.10)
Most m´ar teljes a m´asodrend˝ u v´egesdifferencia m´odszer¨ unk a mindk´et v´eg´en befogott gerenda mozg´asegyenlet´enek megold´as´ara.
5.2.2.
Stabilit´ as
A stabilit´as ebben az esetben is v´egezhet˝o a [18]-ban r´eszletezett Neumann-f´ele m´odszerrel, de a konkr´et kidolgoz´ast ott m´ar nem tal´aljuk meg. Tekints¨ uk a (5.4) s´em´at, de most a jobboldal legyen azonosan nulla. Az i indexeket cser´elj¨ uk le j-re, √ ´es vezess¨ uk be a szok´asos i = −1 jel¨ol´est. V´egezz¨ uk el az ukj = Yjk := λk eiϕk helyettes´ıt´est. Ezek ut´an Yt¯kt,j =
Yjk+1 − 2Yjk + Yjk−1 1 = 2 (λ − 1)2 Yjk 2 τ τ λ
(r´eszletesebben l´asd [18]), valamint k k k k Yj−2 − 4Yj−1 + 6Yjk − 4Yj+1 + Yj+2 1 = 4 Yjk (eiϕ/2 − e−iϕ/2 )4 4 h h 16 k 4 ϕ = 4 Yj sin . h 2
Yx¯kx¯xx,j =
A kapott ´ert´ekeket a s´em´aba behelyettes´ıtve ´es Yjk -val egyszer˝ us´ıtve kapjuk a 1 Eκ2 16 4 ϕ 2 (λ − 1) + sin =0 τ 2λ ρ h4 2 23
n
m
Hiba
Sebess´eg
10
200
0.100539
20
800
0.0274616
3.66107383
40
3200
0.0274616
3.839387076
80
12800
0.00182358
3.922284737
160
51200
0.00046029
3.961806687
5.3. t´abl´azat. A gerenda megval´os´ıt´as´anak hib´aja. karakterisztikus polinomot. Vezess¨ uk be az s :=
sin4 ϕ2 , γ
:=
τ h2
q
Eκ2 ρ
jel¨ol´eseket. A
karakterisztikus polinom gy¨okei ekkor p λ1,2 = 1 − 8γ 2 s4 ± 4γ|s| 4γ 2 s4 − 1. A |λ1,2 | ≤ 1 felt´etel teh´at 4γ 2 ≤ 1 eset´en teljes¨ ul. Ki kell sz˝ urni a t¨obbsz¨or¨os gy¨ok¨ok eset´et, ami csak akkor a´ll fenn, ha 4γ 2 = 1 ´es |s| = 1. Teh´at egy el´egs´eges felt´etel a stabilit´asra 4γ 2 < 1.
5.2.3.
(5.11)
A m´ odszer tesztel´ ese
A kapott m´odszert t¨obb f¨ uggv´ennyel teszteltem, ezek k¨oz¨ ul csak egyet emelek ki: u(x, t) = 12
Eκ2 2 2 t x (1 − x)2 ρ
(x ∈ [0, 1], t ∈ [0, 1])
A pontos megold´ast´ol val´o elt´er´eseket, ´es a konvergenciasebess´egeket mutatja be az 5.3 t´abl´azat. J´ol l´athat´o a m´odszer m´asodrend˝ u konvergenci´aja.
5.3.
A k´ etdimenzi´ os rezon´ ans egyszer˝ u modellj´ enek megval´ os´ıt´ asa
5.3.1.
A v´ egesdifferencia s´ ema fel´ır´ asa
A k´etdimenzi´os esetre a´tt´erve tekints¨ uk a (3.5) mozg´asegyenletet. Legyen most Ω = [a, b] × [a, b] (a, b ∈ R, a < b), amint azt m´ar eddig is feltett¨ uk. Ennek a k¨ovetkez˝o ωh feloszt´as´at fogjuk haszn´alni: ωh := {(xi , yj ) ∈ Ω : xi = a + ih, yj = a + jh, i = 0, 1, . . . , n, j = 0, 1, . . . , n} ahol n ∈ N, h = 1/n. Jel¨olj¨ uk γh -val az ωh ∩ Γ halmazt, azaz a peremen nyugv´o pontokat. Az (xi , yj ) : i, j = 1, n − 1 halmazt a tartom´any bels˝o perem´enek fogom a 24
tov´abbiakban nevezni. Ezen t´ ul u(xi , yj ), (ahol (xi , yj ) ∈ ωh ) r¨ovid´ıtett jel¨ol´ese ui,j . A r¨ovid´ıtett jel¨ol´es ebben az esetben is jel¨olheti ´ertelemszer˝ uen az approxim´aci´ot is. A [0, T ] (0 < T ) id˝ointervallum feloszt´asa a szok´asos. A mozg´asegyenletb˝ol kit˝ unik, hogy a v´egesdifferenci´ak fel´ır´as´ahoz sz¨ uks´eg van egy u ´j approxim´aci´os s´em´ara, mely a ∆2 oper´atorban szerepl˝o vegyes deriv´altat k¨ozel´ıti egy adott (xi , yj ) (r´acs)pontban: 1 · h4 (ui−1,j−1 − 2ui,j−1 + ui+1,j−1 ux¯x¯yy,i,j,h :=
(5.12)
−2ui−1,j + 4ui,j − 2ui+1,j +ui−1,j+1 − 2ui,j+1 + ui+1,j+1 ) A fel´ır´as t¨ordel´es´evel azt pr´ob´altam ´erz´ekeltetni, hogy az x ir´any´ u m´asodrend˝ u approxim´aci´okra alkalmaztam az y ir´any´ ut. Ebb˝ol k¨onnyen l´athat´o, hogy a kapott approxim´aci´o is m´asodrend˝ u, ha u ∈ C 6,6 ([a, b]), mert v´eg¨ ul O(h4 ) hibatagot osztunk h2 -tel. A k¨onnyebb olvashat´os´ag kedv´e´ert az id˝or´eteg megjel¨ol´es´et elhagytam, h is elhagyhat´o. Ezek ut´an a diszkretiz´alt ∆2 oper´ator a k. id˝or´etegben ´ıgy n´ez ki: (∆2 )ki,j,h := ukx¯x¯xx,i,j + 2ukx¯x¯yy,i,j + uky¯yy¯y,i,j Most m´ar minden eszk¨oz¨ unk meg van a diszkretiz´alt egyenlet fel´ır´as´ahoz: ρsukt¯t,i,j +
Es3 k (∆2 )ki,j = fi,j 12(1 − ν 2 )
(5.13)
A kapott diszkretiz´alt s´ema explicit, hisz uk+1 es az im´ent bevezetett i,j -et kifejezve ´ jel¨ol´est haszn´alva kapjuk a m´odszer egy l´ep´es´et: uk+1 i,j :=
τ2 k k−1 f − τ 2 D(∆2 )ki,j + 2uki,j − ui,j ρs i,j
(i, j = 1, 2, . . . , n − 2, k = 1, 2, . . . , m − 1). A k = 1 id˝or´eteget approxim´al´o s´em´at a hull´amegyenlettel ´es az egydimenzi´os esettel anal´og m´odon lehet levezetni ´es hasonl´oan l´athat´o be a m´asodrend˝ us´eg is a diszkr´et ∆2h oper´ator m´asodrend˝ us´eg´et is felhaszn´alva. A nevezett s´ema: u1i,j := u0i,j + τ
∂u τ2 k τ2 (xi , yj , 0) + fi,j − D(∆2 )ki,j ∂t 2ρs 2
(i, j = 1, 2, . . . , n − 2, k = 1, 2, . . . , m − 1). Mivel a ∆2h diszkr´et oper´ator s´em´aja egy ir´anyba 5 pontos, az egydimenzi´os esethez hasonl´oan itt is id˝or´etegen bel¨ uli approxim´aci´ora szorulunk a peremhez legk¨ozelebbi bels˝o pontok kisz´am´ıt´as´ahoz. Erre h´aromf´ele m´odszert is kipr´ob´altam. 25
Els˝ok´ent az egydimenzi´os esetben megismert (5.10) approxim´aci´ot alkalmaztam mindig a norm´alis (azaz a peremre mer˝oleges) ir´anyba. Mivel a sarkok eset´eben k´et ilyen ir´any is lehets´eges, ez´ert ott a kett˝o a´tlag´at vettem, vagy ha u ´gy tetszik, a k´et s´em´at o¨sszeadtam ´es kifejeztem a k´erd´eses pontot. M´asodikk´ent is az (5.10) s´em´at vettem alapul, de most nem csak a norm´alis ir´anyba, hanem a perem vonal´ahoz k´epest 45◦ -os ir´anyba u ´gy, hogy a k´erd´eses pont rajta legyen. Ez ut´obbi lehet˝os´egb˝ol kett˝o is van. ´Igy a mind¨osszesen h´arom approxim´aci´o a´tlag´at vettem. A sarkok eset´eben is h´arom approxim´aci´o ´atlag´at vettem: a k´et norm´alis ir´any´ ut, ´es az a´tl´ora illeszked˝ot. Harmadikk´ent egy u ´j, de m´eg mindig egydimenzi´os approxim´aci´ot vezettem le n´egy pontra. A levezet´es idej´ere t´erj¨ unk vissza az egydimenzi´os jel¨ol´esekhez. A bal oldali perem eset´et vizsg´alva a n´egy pont u0 , u1 , u2 , u3 . Az ut´obbi h´arom pont Taylorsorbafejt´es´et f¨olhaszn´alva ad´odik, hogy 1 1 ∂u u1 = (11u0 + 9u2 − 2u3 ) + h (0) + O(h4 ). 18 3 ∂x Itt u ´jra hivatkozhatunk arra, hogy a peremfelt´etel miatt a jobboldal m´asodik tagja 0 ´es a k¨ovetkez˝o negyedrend˝ u approxim´aci´ot haszn´aljuk: 1 u1 := (11u0 + 9u2 − 2u3 ). 18 Ezt a s´em´at az els˝o approxim´al´ashoz hasonl´o m´odon alkalmaztam, csak a sarkokn´al az ´atl´o ir´any´aban. Lehets´eges m´eg az approxim´aci´ot elhagyni, ´es a peremfelt´eteleket u ´gy kifejezni, hogy a peremen is ´es a ,,bels˝o peremen” is 0 az ´ert´ek. Ez legyen a ,,nulladik” approxim´aci´o.
5.3.2.
Stabilit´ as
Most is vizsg´alhatjuk Neumann-m´odszerrel a stabilit´ast. Tekints¨ uk az (5.13) diszkretiz´aci´ot, megint csak azonosan nulla jobboldallal. A r´acsf¨ uggv´enynek most k´etdimenzi´osnak kell lennie, azaz ukj,l -ba most Yj,lk := λk eiϕ(p)j eiϕ(q)j -t fogjuk behelyettes´ıteni p, q ∈ N, ϕ(r) = 2πrh (r ∈ N) [18]. p az x, q az y tengelyhez tartozik. Ekkor 1 Yt¯kt,j,l = 2 (λ − 1)2 Yj,lk λτ 16 ϕ(p) Yx¯kx¯xx,j,l = 4 Yj,lk sin4 h 2 ϕ(q) 16 Yy¯kyy¯y,j,l = 4 Yj,lk sin4 h 2 1 Yx¯kx¯yy,i,l = 4 Yj,lk (eiϕ(p) − 2 + e−iϕ(p) )(eiϕ(q) − 2 + e−iϕ(q) ) = h 16 k ϕ(p) 2 ϕ(q) Yj,l sin2 sin 4 h 2 2 26
Vezess¨ uk be az s(r) = sin tes´ıt´eseket ´es a
Yj,lk -val
ϕ 2
(r ∈ N), γ =
q
D τ ρs h2
jel¨ol´eseket. Elv´egezve a behelyet-
val´o egyszer˝ us´ıt´est, kapjuk az
ρs
1 16 (λ − 1)2 + D 4 (s2 (p) + s2 (q))2 = 0 2 λτ h
egyenl˝os´eget. Ebb˝ol ad´odik a karakterisztikus polinom egy szebb alakja: 2 2 2 2 2 λ − 2 − 16γ (s (p) + s (q)) λ + 1 = 0 Ennek gy¨okei λ1,2 = 1 − 8γ 2 S 2 ± 4γS ahol S = s2 (p) + s2 (q). A diszkrimin´ans γ ≤
p 4γ 2 S 2 − 1
1 4
eset´en lesz nem pozit´ıv, hiszen S 2
0 ´es 4 k¨oz¨ott minden ´ert´eket f¨olvesz, a 4-et is. Ilyen esetben |λ1,2 | ≤ 1. T¨obbsz¨or¨os gy¨ok γ = 1/4, S = 2 esetben van. Ezeknek f´eny´eben a γ<
1 4
felt´etel el´egs´eges a m´odszer L2 (ωh ) stabilit´as´ahoz.
5.3.3.
A m´ odszer tesztel´ ese
A m´odszert t¨obb f¨ uggv´ennyel teszteltem, ezek k¨oz¨ ul kett˝ore k¨ozl¨ok eredm´enyeket. Az els˝o id˝ot˝ol f¨ uggetlen: u(x, y, t) = (cos(x) − 1.0)(cos(y) − 1.0) (x, y ∈ [0, 2π], t ∈ [0, 2π])
(5.14)
A m´asodik az id˝ot˝ol is f¨ ugg: u(x, y, t) = cos(t)(cos(x) − 1.0)(cos(y) − 1.0) (x, y ∈ [0, 2π], t ∈ [0, 2π])
(5.15)
Alapesetben m id˝or´eteget sz´amoltattam, ha ett˝ol valami´ert elt´ertem, azt k¨ ul¨on jelzem. Az els˝o teszttel, amit k¨ozl¨ok, l´enyeg´eben a ∆2 oper´ator approxim´aci´oj´anak m´asodrend˝ us´eg´et vizsg´alom: Tekints¨ uk az (5.14) tesztf¨ uggv´enyt. Haszn´aljuk a ,,nulladik” approxim´aci´ot. A tartom´anyt az a := −h, b := 2π + h param´eterekkel adjuk meg. Teh´at a tartom´anyt kib˝ov´ıtett¨ uk egy t´err´eteggel. Ekkor a tesztf¨ uggv´eny a bels˝o peremen teljes´ıti a peremfelt´eteleket, a peremen m´ar az ´ert´ek is, ´es az els˝o deriv´alt is pozit´ıv. ´Igy peremfelt´etelk´ent az ott felvett ´ert´ekeket adjuk meg, a deriv´altak most nem j´atszanak szerepet. Azaz a k´et sz´els˝o t´err´etegen megadtuk a pontos ´ert´ekeket. Az n = 12 eset annak felel meg, mintha az n = 10 esetet vizsg´aln´ank, csak a speci´alis peremfelt´etel¨ unkkel. Ez´ert ezzel a feloszt´assal kezdtem, ´ıgy ¨osszem´erhet˝o lesz a t¨obbi eset tesztj´enek kezdeti feloszt´as´aval. Az eredm´enyeket az 5.4 t´abl´azat tartalmazza. 27
n
m
Hiba
Sebess´eg
12
576
0.605734
24
2304
0.113378
5.34260615
48
9216
0.0248762
4.55768968
96
36864
0.0058459
4.25532424
192
147456
0.00141777
4.12330632
5.4. t´abl´azat. A diszkr´et (∆2 )h oper´ator hib´aja. n
m
Hiba
Sebess´eg
12
576
0.605734
24
2304
0.113378
5.34260615
48
9216
0.0248762
4.55768968
96
36864
0.0058459
4.25532424
192
147456
0.00141777
4.12330632
5.5. t´abl´azat. A kib˝ov´ıtett r´eteges m´odszer hib´aja. Ha a m´odszert u ´gy m´odos´ıtjuk, hogy a peremen lev˝o ´ert´ekeket null´azzuk, akkor m´ar egy a val´os´agban is haszn´alhat´o, de els˝orend˝ u m´odszert kapunk. Ezt az esetet az (5.14) f¨ uggv´ennyel teszteltem. Az eredm´enyek az 5.5 t´abl´azatban tal´alhat´oak. Els˝orend˝ u m´odszert kapunk akkor is, ha a nulladik approxim´aci´ot alkalmazzuk a term´eszetes tartom´anyra. Azaz legyen mostant´ol a := 0, b := 2π. Az ekkor kelet´ t˝ kez˝o eredm´enyeket az 5.6 t´abl´azattal mutatom be. Ugy unik, mintha hossz´ ut´avon a mostani eset gyorsabban konverg´alna. M´egis, ha a k´et t´abl´azat els˝o sor´at o¨sszevetj¨ uk, azt l´atjuk, hogy az el˝oz˝o eset hib´aja kisebb. Kipr´ob´altam m´eg egy az el˝oz˝o esettel o¨sszevethet˝o sz´am´ıt´ast a k¨ovetkez˝o param´eterekkel: n = 94, m = 35344. Ekkor 0.224391 hib´at kaptam. Teh´at ebben az esetben is az el˝oz˝o m´odszer gy˝oz¨ott, de u ´gy t˝ unik, az elt´er´es nem sz´amottev˝o. n
m
Hiba
Sebess´eg
10
400
2.13313
20
1600
1.15431
1.84796978
40
6400
0.56445
2.04501727
80
25600
0.2668
2.11562969
160
102400
0.131591
2.02749428
320
409600
0.0654146
2.0116
5.6. t´abl´azat. Az approxim´aci´ot nem alkalmaz´o v´altozat hib´aja.
28
n
m
Hiba
Sebess´eg
10
400
0.26031
20
1600
0.0945783
2.7523
40
6400
0.00276223
3.4240
80
25600
0.00739048
3.7376
160
102400
0.0019137
3.869
320
409600
0.000486841
3.9309
5.7. t´abl´azat. Az els˝o approxim´aci´ot alkalmaz´o v´altozat hib´aja az (5.14) tesztf¨ uggv´enyen. n
m
Hiba
Sebess´eg
10
400
0.173242
20
1600
0.0922018
1.8789
40
6400
0.0280774
3.2838
80
25600
0.00763911
3.6755
160
102400
0.00199065
3.8375
320
409600
0.000508023
3.9184
5.8. t´abl´azat. Az els˝o approxim´aci´ot alkalmaz´o v´altozat hib´aja az (5.15) tesztf¨ uggv´enyen. Ezek ut´an vizsg´aljuk azokat a v´altozatokat, ahol alkalmaztunk approxim´aci´ot a bels˝o peremre. Az els˝o approxim´aci´ot alkalmaz´o v´altozat hib´aj´at mindk´et tesztf¨ uggv´eny eset´en k¨ozl¨om az 5.7, 5.8 t´abl´azatokban. Az (5.15) f¨ uggv´eny eset´en kipr´ob´altam azt is, hogy mi t¨ort´enik, ha az id˝or´etegeket tov´abb sz´amolom. Azt tapasztaltam, hogy a hib´ak egy id˝o ut´an nem n˝onek. A k´ıs´erletet u ´gy v´egeztem, hogy az eddigi t´abl´azatok egyes sorainak megfelel˝o param´eterez´esenk´ent sz´amoltam el˝osz¨or m r´etegig, majd 2m r´etegig, majd 4m r´etegig etc. Addig v´egeztem ezt, am´ıg k´et r´ak¨ovetkez˝o esetben azonos hib´at kaptam. Az utolj´ara kapott, ,,v´egs˝o” hib´akat az 5.9 t´abl´azat tartalmazza. Mindh´arom t´abl´azat azt mutatja, hogy m´odszer¨ unk m´asodrend˝ u, s˝ot nem v´eges id˝o eset´en is tartja ezt. A m´asodik approxim´aci´ot haszn´al´o elj´ar´as hib´aj´at csak az (5.15) f¨ uggv´ennyel t¨ort´en˝o teszt eset´eben k¨ozl¨om az 5.10 t´abl´azatban. A m´odszer hib´aja nem t´ ul je´ lent˝osen, de cs¨okkent. Erdekes, hogy a konvergenciasebess´egre ugyanez mondhat´o el. A harmadik approxim´aci´ot alkalmaz´o v´altozatra vonatkoz´o eredm´enyeket mindk´et tesztf¨ uggv´enyre bemutatom az 5.11 ´es 5.12 t´abl´azatokban. Ezeket o¨sszevetve az 5.7 ´es 5.8 t´abl´azatokkal, meglep˝o eredm´enyt tal´alunk: A magasabbrend˝ u appro29
n
m
Hiba
Sebess´eg
10
400
2.8445
20
1600
0.477454
5.9576
40
6400
0.128435
3.7175
80
25600
0.0343368
3.7404
160
102400
0.000891134
3.8532
5.9. t´abl´azat. Az els˝o approxim´aci´ot alkalmaz´o v´altozat ,,v´egs˝o” hib´aja az (5.15) tesztf¨ uggv´enyen. n
m
Hiba
Sebess´eg
10
400
0.173242
20
1600
0.0922018
1.8789
40
6400
0.0280774
3.2838
80
25600
0.00763911
3.6755
160
102400
0.00199065
3.8375
320
409600
0.000508023
3.9184
5.10. t´abl´azat. Az m´asodik approxim´aci´ot alkalmaz´o v´altozat hib´aja az (5.15) tesztf¨ uggv´enyen. xim´aci´ot alkalmazva nem kapunk jobb eredm´enyt, igaz a konvergenciasebess´eg javul.
n
m
Hiba
Sebess´eg
10
400
0.764936
20
1600
0.168977
4.5269
40
6400
0.0367815
4.5941
80
25600
0.00849313
4.3625
160
102400
0.00205143
4.11
320
409600
0.000504048
4.0699
5.11. t´abl´azat. A harmadik approxim´aci´ot alkalmaz´o v´altozat hib´aja az (5.14) tesztf¨ uggv´enyen.
30
n
m
Hiba
Sebess´eg
10
400
0.85577
20
1600
0.186419
4.5906
40
6400
0.039388
4.7329
80
25600
0.00902204
4.3658
160
102400
0.00216329
4.1705
320
409600
0.000529312
4.087
5.12. t´abl´azat. A harmadik approxim´aci´ot alkalmaz´o v´altozat hib´aja az (5.15) tesztf¨ uggv´enyen.
5.4.
A biharmonikus oper´ ator saj´ at´ ert´ ek probl´ em´ aj´ anak numerikus megold´ asa befogott v´ ekony izotr´ op lemez peremfelt´ etelei mellett
A kor´abbiakban bemutattam a ∆2 oper´atort approxim´al´o diszkr´et ∆2h oper´atort. Az a´ltalunk vizsg´alt Ω tartom´any legyen tov´abbra is a 3.2.4-ben defini´alt n´egyzetes tartom´any, feloszt´asa is legyen az eddigi ω. ´Irjuk f¨ol ∆2 -t m´atrix alakban, azaz h
k´esz´ıts¨ unk egy olyan A m´atrixot, amit ha egy olyan ~b vektorral szorzunk, amiben az u : Ω → R ui,j ´ert´ekei vannak, akkor az eredm´enyben a ∆2h (ui,j ) ´ert´ekek lesznek a megfelel˝o helyeken (u ∈ D4 (Ω), xi , yj ∈ ω, i, j = 0, 1, . . . , n). Ezzel visszavezett¨ uk az oper´ator saj´at´ert´ekprobl´em´aj´at egy m´atrix saj´at´ert´ekprobl´em´aj´ara. V´alasszuk most azt az esetet, amikor az (ui,j )ni,j=0 m´atrix, sorfolytonosan van t´arolva a szorzand´o b vektorban. Ebben a pontban mostant´ol ∆2h jel¨olheti a reprezent´al´o A m´atrixot is. Az eml´ıtett m´atrixba be kell ´ep´ıteni a peremfelt´eteleket is. V´alasszuk a legegyszer˝ ubb esetet, mikor is a nulladik approxim´aci´ot alkalmazzuk. Ekkor a perempontoknak megfelel˝o sorok ´es oszlopok null´ak lesznek. Ha ezef2 -val, a vektorokat a ket elhagyjuk szimmetrikus m´atrixot kapunk, amit jel¨olj¨ unk ∆ h
szok´asos jobbir´any´ u ny´ıllal jel¨ol¨om a v´altoz´ok f¨ol¨ott. A kapott m´atrixokkal fel´ırhat´o a nulladik approxim´aci´ot alkalmaz´o m´odszer¨ unk. Az els˝o id˝or´eteg kisz´am´ıt´asa a k¨ovetkez˝ok´eppen alakul: ~u1 := ~u 0 + τ~u 01 +
τ 2 ~0 τ2 f − D ∆2h~u 0 2ρs 2
(5.16)
A k. id˝or´eteg kisz´am´ıt´asa pedig ´ıgy n´ez ki k = 2, 3, . . . , m. ~u k :=
τ 2 ~k f + (2I − Dτ 2 ∆2h )~u k−1 − ~u k−2 ρs
(5.17)
Mindk´et k´eplet eset´en a vektorok (n + 1)2 , a m´atrixok (n + 1)2 × (n + 1)2 m´eret˝ uek. ~u l -ben az l. id˝or´eteg adott vagy sz´amolt ´ert´ekei vannak (l = 0, 1, . . . , n). ~u 01 a 31
kezdeti ´ert´ekk´ent adott els˝o deriv´altakat, f~ l a jobboldal l. id˝or´eteg eset´en adott f2 ´ert´ekeit tartalmazza sorfolytonosan. I az egys´egm´atrixot jel¨oli. A szimmetrikus ∆ h
2
m´atrix eset´en a k´epletek ugyan´ıgy n´eznek ki, csak a vektorok (n − 3) , a m´atrixok (n − 3)2 × (n − 3)2 m´eret˝ uek. Ebben az esetben az ~u l eredm´enyt term´eszetesen ki kell b˝ov´ıteni a k´et null´akat tartalmaz´o t´err´eteggel. Ezzel a m´odszerrel a v´artaknak megfelel˝oen a m´ar ismertetett nulladik approxim´aci´os esettel azonos eredm´enyeket kaptam. Csak az els˝o tesztf¨ uggv´enyt vizsg´altam, felt´etelezhet˝o, hogy a m´asodikra is azonos eredm´enyeket kapn´ank. Ennek a m´odszernek a t´arig´enye az O(n4 ) m´eret˝ u m´atrix miatt j´oval nagyobb a m´atrix n´elk¨ uli megold´ashoz k´epest (amikor csak O(n2 ) a t´arig´eny, ha nem t´aroljuk az eredm´enyt). Emiatt a gyakorlati jelent˝os´ege kev´es, viszont arra j´o volt, hogy leteszteltem a gener´alt ∆2h m´atrix helyess´eg´et. Ugyanis a k´es˝obbiekben ezt a m´atrixot fogjuk m´eg haszn´alni. L´attuk, hogy a ∆2 biharmonikus oper´ator saj´at´ert´ekprobl´em´aj´at numerikusan megoldhatjuk u ´gy, hogy a ∆2h m´atrix saj´at´ert´ekprobl´em´aj´at oldjuk meg, s˝ot, val´oj´af2 m´atrixot haszn´alnunk. Mivel ez ut´obbi m´atrix szimmetrikus, sokkal ban el´eg a ∆ h
pontosabb numerikus m´odszerek alkalmazhat´oak. Az eml´ıtett m´atrixszal kapott els˝o nyolc saj´atf¨ uggv´enyt n = 20 eset´en az 5.1 a´br´an mutatom be. Az alakok a [12]-ben v´azoltakkal ¨osszhangban a´llnak. A kapott, ´es [12]-ben k¨ozel´ıtett relat´ıv saj´atfrekf2 -val sz´amolt ´ert´ekek venci´akat az 5.13 t´abl´azat tartalmazza. Itt is l´atszik, hogy a ∆ h
´ j´ol k¨ozel´ıtik a m´as m´odszerrel sz´amoltakat. Erdekes megjegyezni, hogy n´emely n param´eterre a m´asodik ´es harmadik saj´atf¨ uggv´eny nem az elv´art alakot veszi f¨ol, hanem olyat, aminek nyugv´o pontjai valamelyik a´tl´ora illeszkednek. A frekvenci´ak viszont ekkor is a helyes ir´anyba v´altoznak. Azt sem tudom meg´allap´ıtani, hogy a feloszt´as s˝ ur˝ us´eg´et˝ol f¨ ugg, mert a probl´ema p´eld´aul n = 10 ´es n = 40 eset´en is el˝oj¨on.
5.5.
A k¨ ozvetlen ¨ osszekapcsol´ as megval´ os´ıt´ asa
A modell fel´ır´asakor haszn´alt jel¨ol´eseket alkalmazom itt is.
5.5.1.
Az egydimenzi´ os eset
A (4.2) egyenletb˝ol diszkretiz´al´assal kapjuk a Eκ2 k Θ ux¯r xr x¯r xr ,r,i = − ukx¯s x¯s ,st,L ρ µ
(k = 1, 2, . . . , m)
(5.18)
s´em´at, ahol ukx¯s x¯s ,st,L = ukx¯s xs ,st,L−1 32
(k = 1, 2, . . . , m)
(5.19)
f2 m´atrix saj´atf¨ uggv´enyei (n=20). A hozz´atartoz´o relat´ıv 5.1. a´bra. A ∆ h saj´atfrekvenci´ak rendre 1.0, 2.0319, 2.0319, 3.0038, 3.6074, 3.6223, 4.5515, 4.5515.
q
n = 10
n = 20
n = 40
[12]-ben k¨oz¨olt ´ert´ek
1.0
1.0
1.0
1.0
2.0
2.0319
2.0378
2.04
λ3 λ1
2.0
2.0319
2.0378
2.04
λ4
2.9748
3.0038
3.0067
3.01
λ5 λ1
3.408
3.6074
3.6453
3.66
λ6
3.4179
3.6223
3.6620
3.67
4.3843
4.5515
4.578
4.58
4.3843
4.5515
4.578
4.58
λ1
q λ1 λ2
q λ1 q
q λ1 q
q λ1 λ7
q λ1 λ8 λ1
f2 m´atrix els˝o nyolc relat´ıv saj´atfrekvenci´aja k¨ 5.13. t´abl´azat. A ∆ ul¨onb¨oz˝o n pah ram´eterekkel, ´es a [12]-ben k¨oz¨olt k¨ozel´ıt˝o ´ert´ekek a befogott n´egyzet alak´ u lemez eset´en.
33
a [17, 18]-ben haszn´alatos jel¨ol´es. Ez a jobboldalon lev˝o approxim´aci´o csak els˝orendben k¨ozel´ıti a h´ ur v´eg´en lev˝o m´asodik deriv´altat, de most ennyivel is megel´egsz¨ unk. uk ukL -val, hisz a k´et pont egy ´es ugyanaz ukr,i -t ´es uks,L -t is form´alisan helyettes´ıthetj¨ a felt´etelez´es¨ unk szerint. Ezek ut´an ukL -t kifejezve kapjuk az ukL
1 = 6Eκ2 + ρh2
Θ µ
Eκ2 Θ k k k k k k (−ur,i−2 + 4ur,i−1 + 4ur,i+1 − ur,i+2 ) + (2ust,L−1 − ust,L−2 ) ρh2 µ (k = 1, 2, . . . , m) (5.20)
s´em´at. Ezt a s´em´at id˝or´etegenk´ent a h´ ur ´es a gerenda u ´j id˝or´eteg´enek kisz´amol´as´at k¨ovet˝oen alkalmazzuk.
5.5.2.
Az egydimenzi´ os eset tesztel´ ese
Tekints¨ uk az ust (x, t) = cos t(cos x − 1) (x ∈ [−π, π]) (5.21)
ur (x, t) = cos t(cos x − 1) (x ∈ [0, 2π]) (t ∈ [0, 2π])
tesztf¨ uggv´enyeket. Az eredm´enyeket az 5.21 t´abl´azat tartalmazza. L´athat´o, hogy a konvergenciasebess´eg t´ uls´agosan oszcill´al ahhoz, hogy m´asodrend˝ us´eg´et elhiggy¨ uk. A probl´ema onnan ad´odik, hogy (5.19) csak els˝orend˝ u k¨ozel´ıt´est ad. Rem´enyteljes viszont, hogy ez a s´ema csak a h´ ur pontjaira t´amaszkodik. Az´ert a´ll´ıtom ezt, mert a h´ ur stabilit´as´anak (5.2) felt´etele sokkal enyh´ebb, mint a gerenda (5.11) felt´etele. A k´et felt´etelb˝ol ad´odik, hogy v´alaszthatjuk a h´ ur h-j´at a gerenda h-j´anak n´egyzet´evel ar´anyosan. A tov´abbi jel¨ol´esbeli u ¨tk¨oz´esek elker¨ ul´ese v´egett ezt a param´etert is indexelj¨ uk a megszokott m´odon att´ol f¨ ugg˝oen, hogy a h´ ur vagy a gerenda feloszt´as´aban j´atszik szerepet. (5.20) levezet´es´ehez hasonl´oan nyerhet˝o egy ezt az esetet is kezel˝o s´ema: ukL =
6Eκ2 ρh4r
1 +
Θ µh2st
·
Θ Eκ2 k k k k k k −u + 4u + 4u − u + 2u − u r,i−2 r,i−1 r,i+1 r,i+2 st,L−1 st,L−2 ρh4r µh2st
(k = 1, 2, . . . , m) (5.22) Azt a´ll´ıtottam, hogy ´elhet¨ unk a hst = O(h2r ) v´alaszt´assal, de a gyakorlatban azt tapasztaltam, hogy enn´el kevesebb is el´eg. A hst = hr /2 v´alaszt´assal m´ar sz´ep eredm´enyeket mutatott a sz´amolt konvergencia. V´altoztassunk egy kicsit az (5.21)
34
n
m
Hiba
Sebess´eg
10
200
0.1587164149
20
800
0.01761144795
9.0121161
40
3200
0.005164685888
3.4099746
80
12800
0.001392766061
3.7082221
160
51200
0.0004036602327
3.4503425
320
204800
0.0001078262967
3.7343616
640
819200
2.78252027·10−5
3.875130
1280
3276800
7.55379752·10−6
3.683604
5.14. t´abl´azat. A k¨ozvetlen egydimenzi´os ¨osszekapcsol´as maxim´alis hib´aja ´es konvergenciasebess´ege az (5.21) tesztf¨ uggv´enyekkel. n
m
Hiba
Sebess´eg
10
400
0.85577
20
1600
0.186419
4.5906
40
6400
0.039388
4.7329
80
25600
0.00902204
4.3658
160
102400
0.00216329
4.1705
320
409600
0.000529312
4.087
5.15. t´abl´azat. A k¨ozvetlen egydimenzi´os ¨osszekapcsol´as maxim´alis hib´aja ´es konvergenciasebess´ege az (5.23) tesztf¨ uggv´enyekkel. teszteseten: ust (x, t) = cos t(cos x − 1) (x ∈ [0, π]) ur (x, t) = cos t(cos x − 1) (x ∈ [0, 2π])
(5.23)
(t ∈ [0, 2π]) L´athat´o, hogy csak a h´ ur tartom´anya v´altozott, a diszkretiz´aci´okor, ekkor h a fel´ere cs¨okken, ha n-t nem v´altoztatjuk. Ennek az esetnek a hib´aj´at tartalmazza az 5.23 t´abl´azat. Ez alapj´an hihet˝o a m´asodrend˝ u konvergencia. Annyit mindenk´epp bizony´ıt, hogy a probl´ema orvosolhat´o.
35
5.5.3.
A k´ etdimenzi´ os eset
Ahhoz, hogy a k´etdimenzi´os esetre a s´em´at viszonylag kompakt m´odon tudjam k¨oz¨olni, bevezetek k´et jel¨ol´est: k 4 2 k 2 u )k \ (∆ pl i,j,h := h (∆ )i,j,h − 20upl,i,j,h
u bkx¯st x¯st ,n,h := h2 ukx¯st x¯st ,n,h − ukst,n,h (k = 1, 2, . . . , m) Az indexekb˝ol h itt is elhagyhat´o. A (4.4) egyenletb˝ol diszkretiz´al´assal kapjuk a −D∆2h ukpl,i,j =
Θ k u µ x¯s x¯s ,st,L
(k = 1, 2, . . . , m)
(5.24)
s´em´at. Itt is elmondhat´o, hogy a jobboldalon az approxim´aci´o csak els˝orend˝ u. Ha a k¨oz¨os pontot ebben az esetben is ukL -val jel¨olj¨ uk ´es beszorzunk h4 µ12ρ(1 − ν 2 )-tel, n´emi a´trendez´es ut´an kapjuk a k 2 2 k 2 2 2 2 \ − Es µ(∆ upl )i,j = Θh 12ρ(1 − ν )b ux¯st x¯st ,n + 20Es µ + ch 12ρ(1 − ν ) ukL 2
(k = 1, 2, . . . , m) egyenl˝os´eget. Ezt ukL -ra kifejezve kapjuk az explicit ukL =
1 2 u )k − Θh2 12ρ(1 − ν 2 )b \ (−Es2 µ(∆ ukx¯st x¯st ,n ) pl i,j 2 2 2 20Es µ + ch 12ρ(1 − ν ) (k = 1, 2, . . . , m) (5.25)
s´em´at. Ezt is id˝or´etegenk´ent a h´ ur ´es a lemez kisz´am´ıt´as´at k¨ovet˝oen alkalmazzuk.
5.6.
A k´ etdimenzi´ os eset tesztel´ ese
A m´odszert a k¨ovetkez˝o tesztf¨ uggv´enyekkel pr´ob´altam: upl (x, y, t) = cos t(cos x − 1)(cos y − 1) (x, y, t ∈ [0, 2π]) ust (x, t) = cos t(cos(x + π) − 1)2
(5.26)
(x, t ∈ [0, 2π])
A lemez numerikus modelljei k¨oz¨ ul a harmadik approxim´aci´ot alkalmaz´o eset´en a hib´akat ´es a konvergenciasebess´egeket az 5.26 t´abl´azat tartalmazza. A teszteset vizu´alis szeml´eltet´es´et mutatja az 5.2 a´bra. Az egydimenzi´os esett˝ol elt´er˝oen, most az (5.19) s´ema els˝orend˝ us´ege nem jut t´ ul nagy ´erv´enyre, legal´abbis a m´ert feloszt´asokon. Ennek oka az lehet, hogy most a m´asodrend˝ u s´ema t¨obb pontra t´amaszkodik, mint az egydimenzi´os esetben. 36
n
m
Hiba
Sebess´eg
10
400
0.466173
20
1600
0.135874
3.4309
40
6400
0.0315055
4.3127
80
25600
0.00718945
4.3822
160
102400
0.00173931
4.1335
320
409600
0.000430698
4.0384
5.16. t´abl´azat. A harmadik approxim´aci´ot alkalmaz´o v´altozat hib´aja az (5.26) tesztf¨ uggv´enyen.
5.2. ´abra. A harmadik approxim´aci´ot alkalmaz´o v´altozat ´es hib´aj´anak vizualiz´aci´oja az (5.26) tesztf¨ uggv´enyen a k = 1528 id˝or´eteg eset´en. A k´ek ´es a s¨ot´etz¨old pontok a sz´amolt ´ert´ekeket, a piros ill. vil´agosz¨oldek a pontos ´ert´ekeket jel¨olik. A sz´amolt ´ert´ekekhez tartoz´oak fedik a t¨obbit.
37
5.7.
Az egyszer˝ u h´ıdon kereszt¨ uli ¨ osszekapcsol´ as megval´ os´ıt´ asa
5.7.1.
A s´ ema ismertet´ ese
Most is haszn´aljuk a modell ismertet´esekor haszn´alt jel¨ol´eseket. A (4.6) egyenl˝os´eg diszkretiz´al´as´ahoz csak az els˝orend˝ u helyderiv´alt approxim´aci´oj´ara van sz¨ uks´eg¨ unk. Tegy¨ uk ezt az ukst,n − ukst,n−1 := (5.27) h els˝orend˝ u k¨ozel´ıt´essel. Ekkor a diszkretiz´alt jobboldalra kapjuk a k¨ovetkez˝o ukx¯s t,n
k fr,i := −ΘG(xr,i )ukx¯s t,n
(i = 0, 1, . . . , n).
s´em´at. Val´oj´aban az i ∈ {l ∈ 0, 1, . . . , n : a + lh < v ∧ w < a + lh} indexekre nem sz¨ uks´eges a konkr´et magval´os´ıt´askor kisz´amolni, mert ott a s´ema nulla ´ert´eket vesz f¨ol. Az egyszer˝ us´eg kedv´e´ert tegy¨ uk fel mostant´ol, hogy valamely c, d indexekre v = a + ch, w = a + dh ´es vezess¨ uk be a p := d − c jel¨ol´est. A h´ ur v´egpontj´aban az u ´j kit´er´est meghat´aroz´o (4.7) egyenl˝os´egben szerepl˝o integr´alt most az o¨sszetett trap´ez szab´allyal k¨ozel´ıts¨ uk [19]: ukst,n
1 := p
! d−1 X 1 k 1 u + uk + uk . 2 st,c l=c+1 st,l 2 st,d
(5.28)
Kihaszn´altam, hogy csak a (4.5) esettel foglalkozunk. A k. id˝or´eteg sz´am´ıt´asa most a k¨ovetkez˝ok´eppen alakul: El˝osz¨or kisz´amoljuk a h´ ur ´ert´ekeit, majd abb´ol frk -t, ezut´an a gerenda pontjait, ´es v´eg¨ ul alkalmazzuk az (5.28) s´em´at a h´ ur gerend´aval o¨sszekapcsolt v´egpontj´ara.
5.7.2.
Tesztel´ es
Mivel arra az estre, amikor a h´ ur az x = 0 pontban van a gerend´ahoz kapcsolva, egyszer˝ ubb volt tesztf¨ uggv´enyeket fel´ırni, ez´ert erre a v´altozatra mutatok be eredm´enyeket. Teh´at a 4.2-ben v´azolt feladat m´odosul annyival, hogy (4.7) helyett most Z ust (0, t) =
b
G(xr )ur (xr , t)dxr . a
38
(5.29)
n
p
m
Hiba
Sebess´eg
10
3
200
0.1080820
20
5
800
0.0294497
3.6700544
40
11
3200
0.00758772
3.8812318
80
21
12800
0.00193223
3.9269238
160
41
51200
0.000487399
3.9643701
5.17. t´abl´azat. A m´odos´ıtott feladat megold´asa az (5.30) tesztf¨ uggv´enyen. legyen igaz. A tesztf¨ uggv´enyek legyenek a k¨ovetkez˝ok: Eκ2 2 2 t x (1 − x)2 (x, t ∈ [0, 1]) ρ µEκ2 ust (x, t) = 12 · Θρ 1 2 1 1 4 2 2 2 3 2 2 2 4 (v + vw + w ) − (v + w)(v + w ) + (v + v w + v w + vw + w ) · 3 2 5 2 √ (x, t ∈ [0, 1]) x + Θµt
ur (x, t) = 12
(5.30) K¨onnyen ellen˝orizhet˝o, hogy ezek a f¨ uggv´enyek a m´odos´ıtott feladatot megoldj´ak, ha a h´ ur a gerenda k¨ozep´ehez van k¨otve, ´es v, w erre a k¨oz´eppontra szimmetrikusan helyezkednek el, ez´ert ´ıgy v´egeztem a teszteket. A teszt sor´an a gerenda a pontos jobboldalt kapta, ´ıgy csak a gerenda h´ urra val´o hat´asa t¨ ukr¨oz˝odik. Az ´ıgy sz´amolt hib´akat ´es konvergenciasebess´egeket mutatja az 5.17 t´abl´azat. p az el˝oz˝o pontban lett bevezetve.
5.8.
Az egyszer˝ u h´ıdon kereszt¨ uli ¨ osszekapcsol´ as megval´ os´ıt´ as´ anak v´ altozata
5.8.1.
A s´ ema ismertet´ ese
Haszn´aljuk 5.7 jel¨ol´eseit. Most csak azokkal a feloszt´asokkal foglalkozunk, amik eset´eben p p´aratlan, Jel¨olj¨ uk l-lel a kapcsol´od´asi pontnak a gerenda diszkretiz´aci´oj´aban felvett index´et. Azaz a gerenda xr,l pontja ,,f¨ol¨ott” a k´epzeletbeli h´ıdra kapcsol´odik a h´ ur. Ekkor c = l − p/2 valamint d = l + p/2, ahol az oszt´as a marad´ekos oszt´ast jel¨oli. Teh´at a h´ ur a h´ıd k¨ozep´en kapcsol´odik hozz´a. (4.9)-t diszkretiz´alva
39
kapjuk: d−1 X h
Gi
i=c
k k fst,n + Θukx¯st x¯st ,st,n fr,i − Eκ2 ukx¯r xr x¯r xr ,r,i = ρ µ
(5.31)
(k = 1, 2, . . . , m). Az integr´alt most a Darboux-f´ele k¨ozel´ıt˝o o¨sszegekkel approxim´altam. ukst,n -t kifejezve ad´odik az d−1
ukst,h
X h2 µ Eκ2 µ k = f − · Θ(p − 1)ρ i=c r,i Θ(p − 1)ρh2
ukr,c−2 − 3ukr,c−1 + 3ukr,c − ukr,c+1 − ukr,d−2 + 3ukr,d−1 − 3ukr,d + ukr,d+1
(5.32)
h2 k − fst,n + 2ukst,n−1 − ukst,n−2 Θ (k = 1, 2, . . . , m) s´ema. Teh´at kaptunk egy u ´j s´em´at a gerend´anak a h´ urra t¨ort´en˝o hat´as´anak kifejez´es´ere.
5.8.2.
Tesztel´ es
A teszteset egyezzen meg a 5.7.2-ban ismertetettel. Azaz most is m´odosul a feladat annyival, hogy a h´ ur x = 0 fel˝oli v´eg´et kapcsoljuk a h´ıdhoz. Erre az esetre a s´ema hasonl´oan vezethet˝o le, mint az el˝oz˝o pontban. A sz´amol´asi hib´akat, ´es a konvergenciasebess´egeket az 5.18 t´abl´azat tartalmazza. Az ATLASZ-on m´ert eredm´enyek az utols´o sorban nagyon elromlanak. Lehets´eges, hogy programoz´asi hiba volt, de az is lehet, hogy a kerek´ıt´esi hib´ak torl´odtak f¨ol. A saj´at g´epemen n = 320 ´es n = 640 sorokhoz tartoz´o hib´ak ar´anya m´eg 4.008 volt. Tov´abb nem sz´amoltattam.
40
n
p
m
Hiba
Sebess´eg
20
9
800
0.0310652
40
17
3200
0.00767899
4.045
80
35
12800
0.00189032
4.0622
160
71
51200
0.000468709
4.033
160
71
51200
0.000468711
320
143
204800
0.000116728
4.0239788
640
287
819200
0.0000298277
3.9134093
1280
575
3276800
0.0000426785
0.6988929
5.18. t´abl´azat. A m´odos´ıtott feladat megold´asa az (5.30) tesztf¨ uggv´enyen a v´altoztatott m´odszerrel. A hatodik sort´ol az ATLASZ-on m´ert eredm´enyek vannak.
41
6. fejezet A tov´ abbi kutat´ omunka lehets´ eges ir´ anyai 6.1.
Pontosabb fizikai modellek haszn´ alata
A k¨ ul¨onb¨oz˝o fizikai modellek ismertet´esekor felv´azoltam n´eh´anyat, amit nem val´os´ıtottam meg. Ezek k¨oz¨ ul mindet ´erdemes haszn´alni. M´egis kiemelek p´arat, amely viszonylag kev´es munka befektet´es´evel megval´os´ıthat´o lehet. ´ Erdemes bevinni csillap´ıt´asi tagokat a h´ ur egyenlet´ebe. Mivel a h´ ur csillap´ıt´asa frekvenciaf¨ ugg˝o, ez´ert k´et tagot is be kell vinn¨ unk: egy els˝o ´es egy harmadrend˝ u id˝oderiv´altat [16]. Ilyen tagokat tartalmaz´o egyenletet tal´alunk [16]-ban. A rezon´ans alakj´anak pontosabb modellez´ese is k¨onnyen megoldhat´o. A legkevesebb v´altoztat´ast az ig´enyli az a´ltalam megval´os´ıtott fizikai modellben, ha nem n´egyzet, hanem t´eglalap alak´ u a lemez. Enn´el nem sokkal t¨obb munk´at ig´enyelne a n´egyzet alak´ u lemez bal f¨ols˝o ´es jobb als´o sark´anak lemetsz´ese egy-egy a sz´elekkel 45 fokot bez´ar´o vonallal. A dupla elhal´as jelens´eg´ehez a h´ ur nem csak a kalap´acs u ¨t´es´evel p´arhuzamos ir´any´ u kit´er´es´et, hanem az arra mer˝oleges kit´er´est is modellezni kell. Ezt r´eszletesen kidolgozva [1]-ben tal´aljuk meg. Erre m´ar nem igaz, hogy kev´es munka ´ar´an megval´os´ıthat´o, ez´ert ez egy t´avolabbi c´el.
6.2.
Pontosabb numerikus modellek, tesztek haszn´ alata
A lemez eset´en a ,,bels˝o perem” k¨ozel´ıt´es´et is lehet jav´ıtani p´eld´aul u ´gy, hogy t¨obb ir´any´ u approxim´aci´ot v´egz¨ unk a m´asodrend˝ u s´ema eset´en is.
42
Az o¨sszekapcsol´asok eset´en az (5.19) approxim´aci´o helyett m´asodrend˝ ut alkalmazva jobb eredm´enyeket ´erhet¨ unk el. Az 5.7.2-ben tal´alhat´o teszt jav´ıt´asa u ´gy, hogy a h´ ur gerend´ara val´o hat´as´anak a hib´aj´at is t¨ ukr¨ozze. Ahol integr´alk¨ozel´ıt˝o formula is szerepel, ott pontosabbak haszn´alata. Lehets´eges az explicit m´odszerek p´arhuzamos´ıt´asa is. Egy-egy id˝or´eteg kisz´amol´as´at oszthatjuk el t¨obb sz´alra. Az o¨sszekapcsol´asok kisz´amol´asa szinkroniz´aci´os pont.
43
7. fejezet Felhaszn´ alt technol´ ogi´ ak, er˝ oforr´ asok, ´ es az azokkal kapcsolatos tapasztalatok A GSL (GNU Scientific Library) tudom´anyos C f¨ uggv´enyk¨onyvt´arat haszn´altam az 5.4-ben k¨oz¨olt megold´asi m´odszer megval´os´ıt´askor. n = 160 param´eter eset´en m´ar nem siker¨ ult lefuttatni a sz´am´ıt´ast. Ami meglep˝o, hogy becsl´eseim szerint el´eg mem´oria ´allt rendelkez´esre. Val´oban a hibakeres´eskor azt tapasztaltam, hogy a m´atrixot inicializ´al´o GSL f¨ uggv´eny nem kiv´etelt okozott (amint azt az el´egtelen mennyis´eg˝ u mem´oria eset´en tenn´e, ´es a programom is kezeln´e ezt az esetet), hanem segmentation fault-ot id´ezett el˝o. Ebb˝ol arra tudok k¨ovetkeztetni, hogy ez a GSL-nek rejtett hib´aja. Azt tapasztaltam, hogy 1Gb-n´al nagyobb ter¨ uletet ig´enyl˝o m´atrix eset´en j¨on el˝o a hiba. A diszkr´et biharmonikus oper´ator saj´at´ert´ekeinek, saj´atvektorainak kisz´am´ıt´as´at GNU Octave-val v´egeztem. Azon bel¨ ul is a ritka m´atrixos reprezent´aci´ot alkalmaz´o elj´ar´asokkal. A programokat, mellyel a t¨obbi sz´am´ıt´ast v´egeztem, C++ nyelven implement´altam. A grafikus fel¨ uletet a QT keretrendszerrel val´os´ıtottam meg, melyb˝ol kihaszn´altam az OpenGL be´agyaz´ast is a h´aromdimenzi´os megjelen´ıt´eshez. T¨obb k´ıs´erletet, tesztet v´egeztem az ELTE ATLASZ (hpc) klaszter´en. Azt tapasztaltam, hogy a saj´at g´epemhez k´epest k¨or¨ ulbel¨ ul k´etszer olyan gyorsan futottak le a sz´am´ıt´asok. Az ott haszn´alt pontosabb, 64 bites aritmetik´ab´ol ad´od´o k¨ ul¨onbs´egek ´eszrevehet˝oek voltak, de nem voltak jelent˝osek. A k¨oz¨olt eredm´enyek a saj´at otthoni g´epemen sz´amoltakat t¨ ukr¨ozi, kiv´eve 5.18 t´abl´azatban foglalt ´ert´ekek egy r´esze. A saj´at g´epem f˝obb param´eterei: K´etmagos Intel(R) Core(TM) Duo T2400 1.83 44
GHz-es CPU, 984.4 MiB DDR2 mem´oria, de a jelen pont els˝o bekezd´es´eben t´argyalt hiba keres´eskor hozz´aadtam m´eg egy 512 MB-os modult. Az ATLASZ egy sz´amol´o node-j´anak f˝obb param´eterei[30]: 2 darab n´egymagos Intel(R) Xeon(TM) E5520 CPU (2.27 GHz), 12GiByte RAM. A fejg´ep eset´en a RAM m´erete 75GiByte. A sz´am´ıt´asok az atlaszon az´ert futhattak az ATLASZ-on dupla sebess´eggel, mert j´oval gyorsabb a mem´oriael´er´ese, ´es a gyors´ıt´ot´arak m´erete is nagyobb. Persze a CPU is kicsit magasabb frekvenci´an fut.
45
8. fejezet Eredm´ enyek ¨ osszefoglal´ asa H´ uros hangszerek, de k¨ ul¨on¨osen a zongora modellez´es´ere siker¨ ult t¨obb matematikai modellt is f¨ol´all´ıtani. A gerend´at tartalmaz´o modell eset´eben h´arom ¨osszekapcsol´asi lehet˝os´eget is ismertettem. A vizsg´alatra kiv´alasztott modelleket diszkretiz´altam v´egesdifferencia m´odszerrel. Siker¨ ult el´egs´eges stabilit´asi krit´eriumokat meghat´arozni az egyes alkatr´eszek diszkretiz´aci´oj´ara. A tesztek seg´ıts´eg´evel a val´os konvergenci´at is m´ertem. Ez is bizony´ıt´asa a m´odszerek gyakorlati alkalmazhat´os´ag´anak. Kider¨ ult, hogy lehets´eges m´asodrend˝ u konvergenciasebess´eget el´erni u ´gy, hogy bizonyos r´eszs´em´ak els˝orend˝ uek, ha k¨ ul¨onb¨oz˝o feloszt´asokat alkalmazunk az egyes o¨sszetev˝okre. R´aad´asul ez u ´gy lehets´eges, hogy nem kell n¨ovelni az id˝or´etegek sz´am´at az azonos feloszt´as´ u esethez k´epest. A k´et esetet konkr´et m´er´essel ¨ossze is hasonl´ıtottam, mely a gyakorlat sz´am´ara kedvez˝o eredm´enyt adott.
46
Irodalomjegyz´ ek [1] Bal´azs Bank: Physics-based Sound Synthesis of String Instruments Including Geometric Nonlinearities, Ph.D. thesis, Budapest University of Technology and Economics Department of Measurement and Infomration Systems, 2006. [2] Bal´azs Bank, Heidi-Maria Lehtonen: Perception of longitudinal components in piano string vibrations, Journal of the Acoustical Society of America, 128(3). [3] Bal´azs Bank, L´aszl´o Sujbert: Generation of longitudinal vibrations in piano strings: From physics to sound synthesis, Journal of the Acoustical Society of America, 117(4), (2005), 2268–2278. [4] Eliane B´ecache, Antoine Chaigne, Gregoire Derveaux, Patrick Joly: Numerical simulation of a guitar, Computers and Structures, 83, (2005), 107–126. [5] Richard E. Berg, David G. Stork: The Physics of Sound, Prentice-Hall, 1982, ISBN 0-13-674283-1. [6] J. Berthaut, M.N. Ichchou, L. J´ez´equel: Piano soundboard: structural behavior, numerical and experimental study in the modal range, Applied Acoustics, 64, (2003), 1113–1136. [7] Malcolm J. Crocker: Handbook of acoustics, John Wiley & Sons, 1998, ISBN 0-741-25293-X. [8] Giuseppe Cuzzocoli, Vincenzo Lombardo: A Physical Model of the Classical Guitar, Including Player’s Touch, Computer Music Journal, 23(2), (1999), 52– 69. [9] M. J. Elejabarrietta, C. Santamaria, A. Ezcurra: Air cavity modes in the resonance box of the guitar: the effect of the sound hole, Journal of Sound and Vibration, 252(3), (2002), 584–590.
47
[10] M. J. Elejabarrietta, C. Santamaria, A. Ezcurra: Fluid-structure coupling in the guitar box: numerical and comparative experimental study, Applied Acoustics, 66, (2005), 411–425. [11] K. A. Exley: Tonal properties of the pianoforte in relation to bass bridge mechanical impedance, Journal of Sound and Vibration, 9(3), (1969), 420–437. [12] Neville H. Fletcher, Thomas D. Rossing: The Physics of Musical Instruments, Springer-Verlag New York, Inc., second edition, 1998, ISBN 0-387-98374-0. [13] N. Giordano: Simple model of a piano soundboard, Journal of the Acoustical Society of America, 102(2), (1997), 1159–1168. [14] N. Giordano: Mechanical impedance of a piano soundboard, Journal of the Acoustical Society of America, 103(4), (1998), 2128–2133. [15] N. Giordano: Sound production by a vibrating piano soundboard: Experiment, Journal of the Acoustical Society of America, 104(3), (1998), 1648–1653. [16] N. Giordano, M. Jiang: Physical Modeling of the Piano, EURASIP Journal on Applied Signal Processing, 7, (2004), 926–933. [17] Stoyan Gisbert, Tak´o Galina: Numerikus m´odszerek 2., ELTE-TypoTEX, 1995, ISBN 963-7546-53-7. [18] Stoyan Gisbert, Tak´o Galina: Numerikus m´odszerek 3., ELTE-TypoTEX, 1997, ISBN 963-7546-77-4. [19] Stoyan Gisbert, Tak´o Galina: Numerikus m´odszerek 1., TypoTEX Kiad´o, 2002, ISBN 963-9326-20-8. [20] Takafumi Hikichi, Naotoshi Osaka: Sound timbre interpolation based on physical modelling, Acoust Sci & Tech, 22(2), (2001), 101–111. [21] Hanna J¨a: Audibilty of the timbral effects of inharmonicity in stringed instrument tones, Acoustics Research Letters Online, 2(3), (2001), 79–84. [22] Adrien Mamou-Mani: Prestress effects on the eigenfrequencies of the soundboards: Experimental results on a simplified string instrument, Journal of the Acoustical Society of America, 131(1), (2012), 872–877. [23] Adrien Mamou-Mani, Jo¨el Frelat, Charles Benainou: Numerical simulation of a piano soundboard under downbearing, Journal of the Acoustical Society of America, 123(4), (2008), 2401–2406. 48
[24] Thomas R. Moore, Sarah A. Zietlow: Interferometric studies of a piano soundboard, Journal of the Acoustical Society of America, 119(3), (2006), 1783–1793. [25] Joseph Derek Morrison, Jean-Marie Adrien: MOSAIC: A Framework for Modal Synthesis, Computer Music Journal, 17(1), (1993), 45–56. [26] Philip M. Morse: Vibration and Sound, American Institute of Physics, paperback edition, 1983, ISBN 0-88318-287-4. [27] Isoharu Nishiguchi: Recent research on the acoustics of pianos, Acoust Sci & Tech, 25(6), (2004), 413–418. [28] A. Stulov: Physical modelling of the piano string scale, Applied Acoustics, 69, (2008), 977–984. [29] Hideo Suzuki: Vibration and sound radiation of a piano soundboard, Journal of the Acoustical Society of America, 80(6), (1986), 1573–1582. ´ am Maulis: Hardver Bemutat´asa, http://www.caesar.elte.hu/hpc/atlasz[30] Ad´ hw.html, 2012.
49
