Numerieke Methode voor Foto-Akoestische Beeldvorming: Een analyse van de k-space methode

Numerieke Methode voor

Foto-Akoestische Beeldvorming

Een analyse van de

k-space methode

Door

W

OUTER

S

WIJGMAN

Begeleid door

D

. D

. C. C. S

Korteweg-De Vries Instituut voor Wiskunde

Faculteit der Natuurwetenschapen, Wiskunde en Informatica UNIVERSITEIT VAN AMSTERDAM

A

BSTRACT

Foto-akoestische beeldvorming is een bio-medische techniek om beelden te kunnen maken van het inwendige van een object. Een onderdeel in het onderzoek naar deze techniek is het oplossen van het voorwaartse probleem. Hier wordt vanuit een begintoestand getracht de golven te simuleren die door het object zullen bewegen. Hoe een golf propageert in een medium kan worden beschreven door de akoestische golfvergelijking. In een homogeen medium is het mogelijk een exacte oplossing voor deze golfvergelijking te geven waarmee dus exact kan berekend worden hoe de golf er op elk tijdstip uit ziet.

In een heterogeen medium zal de golf met een numeriek tijdsstapschema bepaald moeten worden. De toolboxK-WAVEvoorMATLABwaarmee accurate golfsimulaties in heterogene media

gemaakt kunnen worden, maakt gebruik van de k-space methode. Gebaseerd op de akoestische golfvergelijking zal deze k-space methode afgeleid worden die in heterogeen medium nauwkeurig is en computationeel efficiënt. De plaatsafgeleides kunnen namelijk accuraat worden benaderd wegens het gebruik van Fouriertransformaties. Hiermee wordt de vereiste van een dicht raster gematigd ten opzichte van gebruikelijke numerieke schema’s. Deze methode zal geanalyseerd worden en fasefouten voor homogeen en heterogeen medium zullen bestudeerd worden.

Titel: Numerieke Methode voor Foto-Akoestische Beeldvorming; Een analyse van de k-space methode

Auteur: Wouter Swijgman, wouterswijgman@student.uva.nl, 10199020 Begeleider: Dhr. Dr. C. C. Stolk

Einddatum: 16 Juli 2015

Korteweg-De Vries Instituut voor Wiskunde Universiteit van Amsterdam

Science Park 904, 1098 XH Amsterdam

I

NHOUDSOPGAVE

2 Voortplanting van golven 5 2.1 De akoestische golfvergelijking . . . 5

2.1.1 Eendimensionale golfvergelijking in homogeen medium . . . 7

2.2 Eendimensionale oplossing . . . 8 2.2.1 Transportvergelijking . . . 8 2.2.2 Golfvergelijking . . . 9 2.2.3 Beginvoorwaarden . . . 9 2.3 Causaliteit . . . 10 2.3.1 Eén dimensie . . . 10 2.3.2 Twee dimensies . . . 12 2.3.3 Drie dimensies . . . 13

2.4 Twee- en driedimensionale oplossing . . . 14

2.5 Golfkenmerken . . . 14

2.6 Het nut van exacte oplossingen . . . 15

3 Fouriertransformatie 17 3.1 Fourierreeksen . . . 17 3.2 Fouriertransformatie . . . 18 3.2.1 Golfvergelijking transformeren . . . 18 3.2.2 Harmonische oplossing . . . 20 4 Numerieke methodes 21 4.1 Tweede orde eindige differenties . . . 22

INHOUDSOPGAVE

4.2.1 Afleiding k-space methode . . . 24 4.2.2 Vervolg k-space methode . . . 27 4.3 Begincondities implementeren . . . 27

5 Analyse van de methodes 29

5.1 Harmonische golven . . . 30 5.2 Numerieke oplossing in homogeen medium . . . 32 5.3 Numerieke oplossing in heterogeen medium . . . 33

6 Numerieke voorbeelden 35

6.1 Rasterpunten per golflengte . . . 35 6.2 Variabele geluidssnelheid . . . 37 7 Conclusie 39 Populaire samenvatting 41 Foto-akoestische beeldvorming . . . 41 Numeriek schema . . . 42 Bibliografie 45 A Appendix A 47

A.1 Eerste orde stelsel vergelijkingen . . . 47 A.2 k-Space oplossing . . . 50

1

I

NLEIDING

In 1880 was Alexander Graham Bell aan het experimenteren met zijn nieuwe uitvinding: de pho-tophone. Om boodschappen over te dragen, maakte deze uitvinding gebruik van licht. Terwijl hij bezig was merkte hij dat materialen geluid uitzenden wanneer deze worden blootgesteld aan korte lichtpulsen. Dit fenomeen wordt foto-akoestiek genoemd en vindt plaats omdat het materiaal de

Figuur 1.1: Graham Bell en zijn photophone warmte van het licht opneemt, daardoor uitzet

en vervolgens weer inkrimpt wegens afkoel-ing. Door deze beweging ontstaan er drukver-schillen die zich zullen voortplanten in de vorm van drukgolven; de geluidsgolven.[19]Omdat de technologie in die tijd nog niet ver genoeg ontwikkeld was om van dit fenomeen gebruik te kunnen maken, kon er pas honderd jaar later echt mee gewerkt worden. Foto-akoestiek werd namelijk onder andere in de bio-medische sector geïntroduceerd in de vorm van foto-akoestische beeldvorming.[13]

Voor fundamenteel, preklinisch en klinisch onderzoek wordt er intensief gebruik gemaakt

van visualisaties van het inwendige van het geobserveerde.[19]Natuurlijk zijn er al verscheidene technieken bekend waarmee je dergelijke visualisaties kan verkrijgen, zoals de CT-scan en de MRI-scan. Deze bestaande technieken maken echter vaak gebruik van relatief dure apparatuur en/of schadelijke stralingen; alternatieven leveren daarentegen weer een ongedetailleerd beeld op.[13]Een toekomstig alternatief die nu volop in de ontwikkeling is, is de foto-akoestische beeld-vorming en maakt gebruik van het principe dat Bell ontdekte. Foto-akoestische beeldbeeld-vorming, ook wel bekend als PAT (Photo-Acoustic Tomography), is een proces waar geen schadelijke straling of dure apparaten aan te pas komen en eventueel uitgevoerd zal kunnen worden door draagbare apparaten.[19]

Toen in het begin van het nieuwe millennium in vivo beelden werden verkregen met deze techniek, is de ontwikkeling in een razend tempo vooruitgegaan. Het instrumentarium, de

1.1. FOTO-AKOESTISCHE BEELDVORMING

algoritmes voor het reconstrueren van beelden, de visualisatie op verschillende schalen en de in vivo toepassingen in medisch onderzoek gingen met sprongen vooruit.[14]Uiteraard ben ik vooral in dit tweede geïnteresseerd, omdat hier verschillende wiskundige technieken in voor komen.

1.1

Foto-akoestische beeldvorming

Bij de foto-akoestische beeldvorming wordt een een stuk weefsel van een object bestraald door een korte laserpuls waardoor foto-akoestiek zal optreden. De cellen in dit weefsel absorberen namelijk de warmte van de laser en zullen hierdoor uitzetten en vervolgens weer inkrimpen. Op deze manier ontstaan er ultrasone geluidsgolven die zich door het object heen verplaatsen. Deze golven worden uiteindelijk opgevangen door detectoren die rondom het object geplaatst zijn. De laserpuls en de ultrasone detectoren maken een goede combinatie; een duidelijk contrast in het weefsel wordt verkregen doordat sommige cellen beter warmte absorberen dan andere cellen, waarna de ultrasone geluidsgolven zorgen voor een zeer hoge resolutie van de visualisatie.[20] Deze techniek wordt vooralsnog alleen op zacht oppervlakkig weefsel uitgevoerd; bij dieper weefsel zal het licht al door eerder weefsel geabsorbeerd worden, de akoestische golven gaan niet gemakkelijk door bot heen en de golven zullen eerder verzwakken.[13]Een voorname toepassing van deze techniek zal voor bijvoorbeeld borstkankeronderzoek zijn: het is makkelijk doordringbaar weefsel en bovendien absorberen zieke cellen meer warmte dan gezonde cellen (en zenden dus grotere golven uit).[14]Voor verschillende doeleinden kan er tevens gevarieerd worden in verschillende types detectoren. Hoge-frequentie detectoren zullen golven met golflengtes van 10-15 MHz opvangen. Vanwege de kleine golflengtes kan er weliswaar een zeer gedetailleerd beeld opgeleverd worden, maar alleen van oppervlakkig weefsel omdat deze golven snel verzwakt kunnen worden in het medium. Om een beeld te kunnen maken van dieper gelegen delen,

Figuur 1.2: Foto-akoestiek bij cellen

zal gebruikt moeten maken van lage-frequentie detectoren die golflengtes van 2-5 MHz op kunnen vangen. Deze grotere golven worden minder snel verzwakt, maar zullen vanwege de grote golflengtes een minder gedetailleerd beeld ople-veren.[15]Voor de impressie is er van bovenstaande proces een (niet-biologisch verantwoorde) schets gemaakt. Als de detectoren de drukgolven opgevangen hebben begint het in-teressante proces pas. We willen namelijk vanuit deze data te weten komen hoe de begintoestand eruit gezien heeft. Dat is wanneer de wiskunde erbij komt kijken. Het wiskundige deel dat bestudeerd moet worden bestaat eigenlijk uit twee hoofdproblemen; het voorwaartse en inverse probleem.

1.2. WISKUNDIGE TECHNIEK

1.2

Wiskundige techniek

We willen bepalen hoe de druk verandert in het object vanaf het tijdstip dat de laserpuls het object bestraalt tot de golven uitgedoofd zijn. Om dit te kunnen bepalen is een model nodig - beschreven door een verzameling parameters - van het natuurkundige verschijnsel dat optreedt, wat in dit geval de golfbeweging binnen het object is. Daarnaast is er een verzameling observeerbare of meetbare data nodig die afhankelijk is van de modelparameters.

Met het voorwaartse (of directe) probleem probeert men de verzameling meetbare data te voorspellen met behulp van een gegeven verzameling modelparameters (de begintoestand). De gegeven parameters worden in het model gezet en op deze manier wordt er bepaald wat er als data uit moet komen. Om te kunnen bepalen hoe de begintoestand eruit zag als je data hebt verkregen, zal de oplossing voor het voorwaartse probleem als hulpstuk gebruikt worden.

Hier komen we namelijk aan bij het inverse probleem. Bij dit probleem wordt er juist getracht de verzameling modelparameters te bepalen vanuit de verzameling gemeten data en eventueel enkele modelparameters, wanneer het probleem niet voorwaarts opgelost kan worden. Dit kan bijvoorbeeld gedaan worden door de voorwaartse methode te itereren, waar steeds de verzameling modelparameters aangepast wordt totdat de voorspelde data overeenkomt met de gemeten data. Een andere manier is om een expliciete inverse functie op te stellen waardoor je gelijk vanuit je data de modelparameters kan bepalen. Voor verschillende situaties zijn verschillende inverse formules bekend.[8]Toch zal hiervoor ook de voorwaartse methode van belang zijn. Vanuit een begintoestand kan data gesimuleerd worden en via de inverse functie zou deze gesimuleerde data weer moeten resulteren op dezelfde begintoestand.

Er zijn verscheidene onderzoeken waar inverse problemen in voor komen. Om een indruk te geven: er is een tijdschrijft genaamd Inverse Problems die al 30 jaar maandelijks een nieuwe editie uitbrengt met allerlei soorten onderzoek over inverse problemen. Op de website van het tijdschriftwww.iopscience.iop.orgis een overzicht te vinden van alle gepubliceerde artikelen waarvan enkele vrij toegankelijk zijn. Zo zijn er bijvoorbeeld meer dan 300 artikelen die te maken hebben met technieken voor foto-akoestische beeldvorming. Inverse problemen komen echter ook terug in andere contexten, zoals enkele titels duidelijk maken: A tutorial on inverse problems for anomalous diffusion processes, Inverse problems in the Bayesian framework, Joint reconstruction of PET-MRI by exploiting structural similarity, Complex eigenvalues and the inverse spectral problem for transmission eigenvalues.

1.3

k-Wave

Opwww.k-wave.orgis een toolbox te vinden waarmee golven met een zekere nauwkeurigheid gesimuleerd kunnen worden. Het is gebaseerd op een numerieke methode die beschreven en bestudeerd wordt in allerlei artikelen die tevens op de website te vinden zijn. Aan de hand van de artikelen van Treeby et al.[1] en Tabei et al.[2] heb ik gekeken welke numerieke methode

1.3. K-WAVE

er gebruikt wordt. Hoe komt namelijk een dergelijke methode tot stand? Er bestaan uiteraard ook andere methodes; waarom is deze methode geschikter dan andere bestaande methodes? We zullen de methode analyseren in een homogeen en heterogeen medium en bekijken hoe vanuit begincondities een golf gesimuleerd kan worden. Omdat het een numerieke benadering is, zullen er waarschijnlijk ook fouten ontstaan. Een belangrijke vraag is dan natuurlijk hoe exact de methode is voor verschillende media en of dit nog valt te verbeteren.

Omdat bij het proces van het voorwaartse probleem al veel vragen en wiskundige technieken opduiken, ben ik niet aan het inverse probleem toegekomen. Om te begrijpen hoe het numerieke schema werkt dat gebruikt wordt bijK-WAVEzal eerst fundamentele kennis over golven en de

wiskundige technieken vergaard moeten worden. Om het kort samen te vatten is ten eerste een model nodig dat beschrijft hoe een golf zich voortplant. Enkele eigenschappen van het gedrag van golven zullen besproken worden om meer over deze golven te weten te komen. Vervolgens zal ik werken naar de numerieke methode om een golf te kunnen simuleren vanuit een begintoestand en bekijken hoe nauwkeurig deze is. In een medium waar de geluidssnelheid constant is, kan er ook een exacte oplossing geconstrueerd worden waarmee de numerieke simulatie vergeleken kan worden. In het begin zal de scriptie een natuurkundig karakter hebben waar de modellen en eigenschappen op fysische wijze beargumenteerd worden. Toch zal verderop het wiskundige karakter terugkomen wanneer de numerieke methode bestudeerd wordt.

2

V

OORTPLANTING VAN GOLVEN

We zullen deze scriptie beginnen met de fundamentele wiskunde, of eigenlijk natuurkunde, waarop de gehele verdere theorie op berust. We beginnen namelijk door te kijken naar een model dat beschrijft hoe akoestische golven propageren in een medium. Omdat dit een belangrijk model is, zal de totstandkoming ervan bekeken worden vanuit natuurkundige principes. Het propageren van akoestische golven kunnen we beschrijven met een tweede orde partiële differentiaalver-gelijking en met een stelsel van eerste orde verdifferentiaalver-gelijkingen. Deze twee modellen zullen worden afgeleid met behulp van de educatieve webpagina van Stanford[16]. De vergelijkingen kunnen opgelost worden voor de onbekende functie(s) om te weten te komen hoe de golf eruit ziet en verandert in tijd. Hiervoor zijn voornamelijk de boeken van Olver[5]en Strauss[6]bestudeerd. Voor het eendimensionale geval kan een algemene oplossing gegeven worden. Wanneer er zekere begincondities zijn kunnen er expliciete oplossingen gegeven worden. Verder zal er nog naar het gedrag van de voortplanting van golven gekeken worden, om vervolgens naar oplossingen in het meerdimensionale geval te kunnen redeneren. Laten we eerst de afleiding van de twee verschillende modellen bekijken en de relatie ertussen.

2.1

De akoestische golfvergelijking

Akoestische golven zijn longitudinale golven die bestaan uit opvolgingen van drukpulsen of elastische verplaatsingen van een materiaal. Dit betekent dat de uitwijking van de golf in dezelfde richting staat als de verplaatsing ervan. De snelheid van een akoestische golf in een medium is afhankelijk van twee (lokale) eigenschappen van het medium; de dichtheid en de onsamendrukbaarheid. De dichtheid geeft aan hoeveel deeltjes er per volume-eenheid aanwezig zijn en de onsamendrukbaarheid geeft aan in welke mate het medium weerstand beidt tegen samendrukking. Akoestische golven die hoorbaar zijn voor mensen worden geclassificeerd als geluidsgolven; golven met hogere frequentie worden ultrasone geluidsgolven genoemd. Aan de hand van de akoestische golfvergelijking kunnen we dergelijke golven beschrijven:[4]

(2.1) 1 κ(x) ∂2_u ∂t2(x, t) − ∇ · ³ 1 ρ(x)∇u(x, t) ´ = f (x, t).

2.1. DE AKOESTISCHE GOLFVERGELIJKING

Hier duidt u het drukveld,_ρde dichtheid van het medium,_κde onsamendrukbaarheid of Bulk modulus van het medium en f een uitwendige kracht aan. In de vergelijking wordt de gradient -de afgelei-des naar x, y, z - aangeduid met ∇. Divergentie wordt aangeduid met ∇· en geeft aan hoe een vectorveld verandert; dit geeft de netto snelheid van massaverandering van de vloeistof dat van een punt wegstroomt aan per volume-eenheid. Waarom deze verschillende termen de akoestische golf beschrijven zal ik uitleggen door te beginnen met een eendimensionaal geval; deze kan vervolgens uitgebreid worden naar een meerdimensionaal geval, maar de natuurkunde erachter blijft behouden.

Het eerste natuurkundige principe waar de vergelijking op berust is de tweede wet van Newton: F = m · a. Deze wet zegt dat de kracht op een voorwerp gelijk is aan de verandering van de impuls - in deze context wil dat zeggen dat een klein volume in een vloeistof zal accelereren als er een kracht op uitgeoefend wordt.[16]Deze kracht ontstaat door drukverschillen aan de verschillende kanten van het stukje volume. We kunnen de vergelijking massa × acceleratie = kracht = −gradient druk bekijken voor het een-dimensionale geval, waar v de flow van de snelheid van de vloeistof aanduidt:[16]

ρ(x) ·∂v

∂t(x, t) = −

∂u

∂x(x, t).

Het tweede natuurkundige proces dat plaatsvindt is de wet van behoud van massa bij aan-wezigheid van samendrukking en verandering van volume. Als de snelheidsvector v op een bepaalde positie x kleiner is dan op een klein stukje verder - aangeduid met x +∆x - dan noemen we de flow divergerend. Dit betekent dat alle deeltjes binnen dit stuk bewegen, met de deeltjes rond x +∆x sneller dan de deeltjes op x. Hieruit volgt dat het kleine volume uitbreidt, maar hierdoor bevinden zich op dat stuk minder deeltjes per volume-eenheid; de druk van dit kleine volume daalt. De grootte van vermindering van de druk staat in verhouding met een eigenschap van de vloeistof, de onsamendrukbaarheid_κ. De wet van behoud van massa geeft in deze con-text de verhouding vermindering druk = onsamendrukbaarheid × divergentie snelheid, wat de volgende vergelijking geeft:[16]

−∂u

∂t(x, t) =κ(x) ·

∂v

∂x(x, t).

Voor het meer-dimensionale geval geldt exact hetzelfde, behalve dat de afgeleide van u nu naar twee of drie richtingen is - de gradient van u. In de algemene vorm van de golfvergelijking wordt ook een uitwendige kracht in de vergelijking gegeven. In deze scriptie wordt deze gelijk aan 0 genomen en zal om die reden zonder uitleg geïntroduceerd worden. De voorgaande twee formules kunnen uitgebreid worden naar meer dimensies inclusief een uitwendige kracht die druk uitoefent op een klein volume. De golfvergelijking kan dan worden beschreven door het stelsel eerste orde vergelijkingen:[2]

1 κ(x) ∂u ∂t(x, t) = −∇ · v(x, t) + F(x, t), (2.2) ρ(x)∂v ∂t(x, t) = −∇u(x, t). (2.3)

2.1. DE AKOESTISCHE GOLFVERGELIJKING

Om uiteindelijk op de akoestische golfvergelijking (2.1) te komen zullen we (2.3) in (2.2) sub-stitueren. Daarvoor zal eerst de vergelijking als volgt omgeschreven worden:

ρ(x)∂v ∂t(x, t) = −∇u(x, t), ∂v ∂t(x, t) = − 1 ρ(x)∇u(x, t), ∇ ·∂v ∂t(x, t) = −∇ · ³ 1 ρ(x)∇u(x, t) ´ .

Als we vervolgens de andere vergelijking in het stelsel ook omschrijven, kunnen we met behulp van substitutie vanuit bovenstaande vergelijkingen de tweede orde akoestische golfvergelijking (2.1) verkrijgen: 1 κ(x) ∂u ∂t(x, t) + ∇ · v(x, t) = F(x, t), 1 κ(x) ∂2_u ∂t2(x, t) + ∇ · ∂v ∂t(x, t) = ∂ ∂tF(x, t), 1 κ(x) ∂2_u ∂t2(x, t) − ∇ · ³ 1 ρ(x)∇u(x, t) ´ = f (x, t).

Golven zijn dus op twee verschillende manieren te beschrijven: als een tweede orde partiële differentiaalvergelijking (2.1) en als een stelsel van eerste orde vergelijkingen (2.2) en (2.3). Ze beschrijven echter precies hetzelfde en dat geeft dus de vrijheid om degene te gebruiken die in de gegeven situatie het prettigst is om mee te werken. In de numerieke methode die later gepresenteerd wordt zal de tweede orde partiële differentiaalvergelijking het uitgangspunt zijn, waardoor we in dit hoofdstuk dieper op deze differentiaalvergelijking in zullen gaan.

2.1.1 Eendimensionale golfvergelijking in homogeen medium

Laten we aannemen dat de onsamendrukbaarheidκen de dichtheidρvan een zeker medium (bijvoorbeeld water) constant zijn en dat er geen uitwendige krachten uitgeoefend worden op dit medium. We bekijken dan een eendimensionale akoestische golfvergelijking in een homogeen medium. Als bovengenoemde condities in (2.1) geïmplementeerd worden, wordt de volgende vergelijking verkregen: 1 κ ∂2_u ∂t2(x, t) − ∂ ∂x 1 ρ ∂u ∂x(x, t) = 0

Nu de onsamendrukbaarheid en de dichtheid niet meer afhankelijk zijn van een variabele, kan de_ρten eerste door de differentiaaloperators gehaald worden en kan er vervolgens met_κworden vermenigvuldigd om op de volgende vergelijking uit te komen:

(2.4) ∂

2_u

∂t2(x, t) − c 2∂2u

2.2. EENDIMENSIONALE OPLOSSING

waar c =qκ_ρ de snelheid van een geluidsgolf door het medium aangeeft[4] – in het geval van water is dat 1500 m/s. Dit is de eendimensionale versie van de golfvergelijking in een vloeistof waarvan de dichtheid en onsamendrukbaarheid constant zijn en is dus meest eenvoudige variant. Van deze vergelijking willen we vervolgens weten hoe de oplossing u er uitziet; dan kunnen we namelijk voor elke positie en elk tijdstip op een zeker gebied bepalen wat de druk zal zijn.

2.2

Eendimensionale oplossing

Als we in één dimensie een medium bekijken waar de onsamendrukbaarheid en dichtheid constant zijn (2.4), kan er een algemene oplossing voor de golfvergelijking bepaald worden. Wanneer er ook nog begincondities van de golf gegeven zijn kan de exacte oplossing bepaald worden, bekend als d’Alemberts oplossing. Omdat er in deze dimensie een verband is tussen de transportvergelijking en de golfvergelijking, bekijken we om te beginnen oplossingen van de eerste.

2.2.1 Transportvergelijking

We beginnen met de aanname dat de geluidssnelheid c in het medium constant is. De eendimen-sionale transportvergelijking modelleert het transport van een substantie in een medium en wordt gegeven door de eerste orde partiële differentiaalvergelijking:[6]

∂u

∂t(x, t) + c

∂u

∂x(x, t) = 0.

Vervolgens kanξ= x − ct geïntroduceerd worden als de karakteristieke variabele[6]; voor een object met positie x, is_ξde positie van het object relatief tot de observeerder die met snelheid c naar rechts beweegt. Op deze wijze kan het coördinatenstelsel van ruimte-tijd (x, t) omgeschreven worden tot het stelsel met de bewegingscoördinaten (ξ, t). Voor een zekere functie v geldt dan u(x, t) = v(x − ct, t) = v(ξ, t). Als we vervolgens v(ξ, t) in de transportvergelijking zetten, houden we wegens de kettingregel alleen het volgende over:[6]

∂v

∂t(ξ, t) = 0.

Hieruit kan worden opgemaakt dat u(x, t) de transportvergelijking oplost dan en slechts dan als v(ξ, t) bovenstaande vergelijking oplost. Als we deze vergelijking nauwkeurig bekijken zien we dat v alleen afhankelijk is vanξ_{. Er volgt dat u(x, t) = v(}ξ_{) = v(x − ct), en dus is de oplossing voor} de transportvergelijking ook alleen afhankelijk van de karakteristieke variabele. De oplossing is dus een naar rechts bewegende golf met snelheid c.

2.2. EENDIMENSIONALE OPLOSSING

2.2.2 Golfvergelijking

Omdat differentiatie een operator is, kunnen we de golfvergelijking (2.4) met constante c op de volgende manier omschrijven:

∂2_u ∂t2(x, t) − c 2∂2u ∂x2(x, t) = ³∂2 ∂t2− c 2 ∂2 ∂x2 ´ u(x, t) =³∂ ∂t+ c ∂ ∂x ´³∂ ∂t− c ∂ ∂x ´ u(x, t),

waar bij de laatste gelijkheid gebruik gemaakt wordt van het feit dat de kruistermen tegen elkaar wegvallen. Duidelijk is dat u(x, t) een oplossing voor de golfvergelijking is als¡∂

∂t− c∂x∂¢u(x, t) = 0.

Deze vergelijking is echter gewoon de transportvergelijking, maar dan met snelheid −c, dus terugwaarts. Zoals hierboven aangetoond, wordt deze opgelost door u(x, t) = p(x + ct) voor p ∈ C2 een tweemaal continu differentieerbare functie; de oplossing is een golf die naar links beweegt met snelheid c.

Aangezien we de twee operatortermen kunnen omwisselen, zien we dat andere vergelijking ¡∂

∂t+ c∂x∂¢u(x, t) = 0 ook een oplossing voor de golfvergelijking geeft. Dit is de gewone

trans-portvergelijking waarvoor we weten dat de oplossing van de vorm u(x, t) = q(x − ct) is, voor een zekere q ∈ C2; een golf die naar rechts beweegt met snelheid c. Wegens lineariteit van de golfver-gelijking is de som van twee oplossingen wederom een oplossing; we zien dus dat we oplossingen kunnen maken bestaande uit links- en rechtsbewegende golven. De volgende stelling zegt dat zelfs elke oplossing voor de golfvergelijking op deze manier gegeven kan worden. Het bewijs ervan maakt wederom gebruik van karakteristieke variabelen en laat ik achterwege.[5]

Stelling 2.1. Elke oplossing voor de eendimensionale golfvergelijking in homogeen medium kan

worden geschreven als de som van links- en rechtsbewegende golven:

(2.5) _{u(x, t) = p(x + ct) + q(x − ct)}

waar p en q C2-functies zijn.

2.2.3 Beginvoorwaarden

Wanneer er ook beginvoorwaarden gegeven zijn kan een explicietere oplossing in termen van de beginvoorwaarden gegeven worden, bekend als d’Alemberts oplossing. In dit geval spreekt men dan van een beginwaardeprobleem, beter bekend als de Engelse versie initial value problem (IVP). Deze beginvoorwaarden bestaan uit de initiële druk en initiële verandering van de druk en worden beschreven door:

u(x, 0) = f (x), ∂u

∂t(x, 0) = g(x).

We mogen volgens bovenstaande stelling aannemen dat onze oplossing van die vorm is. Als we vervolgens voor deze vorm de beginvoorwaarden bekijken, krijgen we:

2.3. CAUSALITEIT

Om deze twee vergelijking op te lossen voor p en q zal eerst f gedifferentieerd worden en g door c gedeeld worden:

f0_{(x) = p}0_{(x) + q}0(x), 1

cg(x) = p

0_{(x) − q}0_(x).

We kunnen p uitdrukken in g en f door de vergelijking bij elkaar op te tellen, waar na intergratie de uitdrukking volgt: 2p0_{(x) = f}0_{(x) +}1 cg(x), p(x) =1 2f (x) + 1 2c Z x 0 g(z) d z + k.

Nu kan ook q uitgedrukt worden op de volgende wijze:

q(x) = f (x) − p(x) =1 2f (x) − 1 2c Z x 0 g(z) d z − k.

Deze twee uitdrukking kunnen we tot slot weer invullen in de algemene oplossing om d’Alemberts oplossing te verkrijgen: u(x, t) = p(x + ct) + q(x − ct) =1 2f (x + ct) + 1 2c Z x+ct 0 g(z) d z + k + 1 2f (x − ct) − 1 2c Z x−ct 0 g(z) d z − k = f (x + ct) + f (x − ct) 2 + 1 2c Z x+ct x−ct g(z) d z.

2.3

Causaliteit

Een belangrijke eigenschap van een golf is de causaliteit. Aan de hand van het boek van Strauss[6] zal gekeken worden wat causaliteit is en hoe dit in verband staat met de oplossing. Met causaliteit wordt gezegd dat een golf slechts een bepaald domein kan beïnvloeden, maar ook dat een punt alleen beïnvloed kan zijn door een bepaald domein. Dit is tevens een handige manier om subtiel de overstap te maken naar de twee- en driedimensionale oplossing vanuit de eendimensionale oplossing. Van de meerdimensionale oplossing zal ik namelijk geen expliciete afleiding geven, maar causaliteit zal wel helpen de formules intuïtief te begrijpen. We zullen wederom kijken naar een homogeen medium, dus met constante snelheid c.

2.3.1 Eén dimensie

Als er geen initiële verandering van de druk aanwezig is, plant een golf zich in beide richtingen met snelheid c op dezelfde manier voort. Stel dat de golf op plaats x0 begint, dan kunnen we

schetsen waar de golffront zich bevindt op verschillende tijdstippen. Stel dat er op (x0, 0) een

beginvoorwaarde is. Deze beginvoorwaarde kan dan alleen de punten x ≤ x0+ ct en x ≥ x0− ct

beïnvloeden omdat de golf met snelheid kleiner of gelijk aan c in beide directies voortplant, zoals we in Figuur 2.1 zien. We kunnen deze voortplanting ook in het (x, t)-vlak bekijken, waar we

2.3. CAUSALITEIT t = 0 x0 t = 1 x0− c x0+ c t = 2 x0− 2c x0+ 2c

Figuur 2.1: Causaliteit in één dimensie

dan een zogeheten domein van invloed kunnen aangegeven als het gearceerde gebied; hierin bevinden zich dus de punten waarvoor geldt |x0− x| ≤ ct. Het gebied wordt ingesloten door de

lijnen |x0− x| = ct; de golven die zich met de maximale snelheid c voortplanten. Hieruit volgt

dan ook dat u(x, t) = 0 voor |x| > M + t als de beginvoorwaarden voor |x| > M verdwijnen. Met

x t

x + ct = x0 x − ct = x0

(x₀, 0)

Figuur 2.2: Domein van invloed voor het punt (x0, 0)

deze denkwijze kunnen we ook het domein van afhankelijkheid bepalen. Neem hiervoor een vast punt (x0, t0) op een bepaalde plek en tijdstip. Van u(x0, t0) weten we dat deze alleen afhankelijk

is van de begincondities, namelijk van de initiële druk op x_{0− ct0} en x_{0+ ct0}en de beginsnelheid op het interval [x0− ct0, x0+ ct0]. Het interval van afhankelijkheid van het punt (x0, t0) is dan

[x0− ct0, x0+ ct0], en het domein van afhankelijkheid is in de figuur gegeven door het gearceerde

gebied. Duidelijk is te zien dat de karakteristieken_{ξ= x − ct en η= x + ct een belangrijke rol} spelen in het domein van invloed en het domein van afhankelijkheid; het zijn namelijk de randen van de domeinen.

x t

(x0− ct0, 0) (x0+ ct0, 0)

(x0, t0)

2.3. CAUSALITEIT

2.3.2 Twee dimensies

We zagen net in het eendimensionale geval dat punten alleen afhankelijk zijn van beginpunten in een bepaald interval, die ingesloten wordt door de karakteristieken. In het twee- en driedimen-sionale geval zijn deze karakteristieken weer net zo belangrijk, maar nu zijn het oppervlakken in plaats van lijnen. Om dit in te zien kunnen we eerst weer bekijken waar de golffront zich bevindt als deze met maximale snelheid c zich voortplant. In het tweedimensionale geval zien we dat de golffront een cirkel is dat zich uitspreidt, met radius ct. Het karakteristieke oppervlak dat hierbij

t = 0 (x₀, y₀) t = 1 2c t = 2 4c

Figuur 2.4: Causaliteit in twee dimensies

hoort kunnen we op de volgende manier construeren. We nemen een karakteristieke lijn in één dimensie door het punt (x0, t0) en draaien deze om de as van x = x0. Voor de tweedimensionale

golf is er dan het driedimensionale karakteristieke oppervlak:

|x − x0| = q (x − x0)2+ (y − y0)2 = c|t − t0|. x t y (x0, y0, t0)

2.3. CAUSALITEIT

Het karakteristieke oppervlak voor de tweedimensionale golf wordt dan het oppervlak van een dubbele kegel, wat de karakteristieke kegel of lichtkegel van (x0, t0) voorstelt. Deze

karakter-istieke kegel bestaat uit de vereniging van de “toekomstige” en “verleden” kegel; het domein van invloed en het domein van afhankelijkheid. Wederom geldt hier weer dat de oplossing u(x0, y0, t0)

afhankelijk is van beginwaardes binnen de kegel,p_{(x − x}0)2+ (y − y0)2 ≤ ct0.

2.3.3 Drie dimensies

Voor de golf in drie dimensies hebben we precies weer hetzelfde principe, alleen nu is de positie van de golffront een sfeer waarvan de radius met snelheid c groter wordt als er geen beginsnelheid aanwezig is. Stel namelijk dat we een punt (x0, y0, z0) hebben waar de golf ontstaat, dan spreidt

deze zich in alle richtingen met de zelfde snelheid c uit. Ook voor deze golfbeweging hebben dan

t = 0 (x0, y0, z0) t = 1 2c t = 2 4c

Figuur 2.6: Causaliteit in drie dimensies

weer een karakteristiek oppervlak, maar dat is in dit geval een vierdimensionale hyperkegel gegeven door:

|x − x0| =

q

(x − x0)2+ (y − y0)2+ (z − z0)2 = c|t − t0|.

Voor de golf in één en twee dimensies zagen we dat de oplossing u(x0, t0) afhankelijk was

beginwaardes op het gebied |x0− x| ≤ ct0. Als we dit uitbreiden naar onze driedimensionale golf,

zouden we verwachten dat de oplossing u(x0, t0) afhangt van de beginwaardes op de punten

in de bol |x0− x| ≤ ct0. Dit is echter niet het geval. Volgens de principe van Huygens geldt in

drie dimensies dat de waarde u(x0, t0) alleen af van de beginwaardes voor x op het oppervlak

van de bol, |x − x0| = ct0, en dus niet van de waardes binnen de bol! Uit deze principe kunnen

we ook opmaken dat beginwaardes op een punt x0 alleen de punten op het oppervlak van de

lichtkegel |x0− x| = ct beïnvloeden. Dit betekent dus dat elke oplossing van de golfvergelijking in

2.4. TWEE- EN DRIEDIMENSIONALE OPLOSSING

2.4

Twee- en driedimensionale oplossing

We zagen net dat de waarde u(x0, t0) van de driedimensionale golfvergelijking alleen afhankelijk

is van beginwaardes op de sfeer S, gedefinieerd door |x0− x| = ct. Stel dat we de volgende

beginwaardes hebben:

u(x, 0) = f (x), ∂u

∂t(x, t0) = g(x).

We zien in Kirchoff ’s formule voor de oplossing van de driedimensionale golfvergelijking weer terug dat deze alleen afhankelijk is van beginwaardes op de sfeer S:[6]

u(x₀, t_{0) =} 1 4πc2_t 0 Ï S g(x) dS + ∂ ∂t0 ³ 1 4πc2_t 0 Ï S f (x) dS´.

Gezien de omvang van de afleiding voor deze formule laat ik die achterwege.[6] Vanuit deze oplossing kunnen we vervolgens naar een oplossing van de tweedimensionale golfvergelijking werken. Met behulp van de method of descent kan vanuit de driedimensionale oplossing naar de tweedimensionale oplossing gewerkt worden. Deze afleiding laat ik echter eveneens achter-wege.[17]De oplossing van de tweedimensionale golfvergelijking met bovenstaande begincondities wordt dan:[6] u(x₀, y0, t0) = 1 2πc Ï D g(x₀, y₀) q c2_t2 0− (x − x0)2− (y − y0)2 dx d y + 1 2πc ∂ ∂t0 Ï D f (x, y) q c2_t2 0− (x − x0)2− (y − y0)2 dx d y

In deze vergelijking wordt geïntergreerd over de schijf D; dit zijn alle punten die voldoen aan de ongelijkheidp_{(x − x}0)2+ (y − y0)2 ≤ ct0 . Deze oplossing is dus afhankelijk van punten binnen

de kegel en niet alleen van punten op de rand, zoals al opgemerkt is bij de causaliteit bij een tweedimensionale golf.

2.5

Golfkenmerken

We zullen dit hoofdstuk afsluiten met nog enkele begrippen die een bij een golf van belang zijn. Als voorbeeld kan de vlakke golf u(x, t) = sin(x − ct) bekeken worden; een sinusgolf die zich naar rechts verplaatst. Als deze vastgezet wordt op bijvoorbeeld een tijdstip t0= 0, kan deze

afgebeeld worden als in Figuur 2.7. De golf kan vanaf een bepaald punt gevolgd worden over het ruimtelijke domein, en zal zich na een bepaalde afstand herhalen; een dergelijke ruimtelijke herhaling is een golfcykel. De faseφvan een periodieke golf is de fractie van de golfcykel die is afgelegd met betrekking tot de oorsprong (waarφ_{= 0). Als deze golf vergeleken wordt met} een golf die bijvoorbeeld met π₂ verschoven is, dan kan het verschil van π₂ gezien worden als een faseverschil van 1₄, ofwel 25%. De golflengte_λis de lengte van de golfcykel; de afstand tussen twee opeenvolgende punten met dezelfde fase. De golf zal zich met snelheid c naar rechts toe

2.6. HET NUT VAN EXACTE OPLOSSINGEN

Figuur 2.7: De golf u(x, t) = sin(x − ct) in het (u, x)-vlak.

verplaatsen en zal na een bepaalde tijd er weer hetzelfde uitzien. Om precies te zijn gebeurt dit op T = λ_c en op deze wijze geeft T de periode van de golf weer. Tot slot zijn er nog twee kenmerken die bij de Fourieranalyse van belang zullen zijn. Het golfgetal k is het aantal radialen per afstandseenheid en wordt gegeven door k =2_λπ. Met dit getal wordt dus aangegeven hoeveel golfcykels er op een bepaalde afstand kunnen worden afgelegd. Tot slot is het frequentiegetal_ω het aantal radialen per tijdseenheid en wordt gegeven doorω₌2_Tπ.

2.6

Het nut van exacte oplossingen

Hoewel in deze scriptie het simuleren van een golf met behulp van numerieke benaderingen voor de golfvergelijking juist het doel is, is het toch ook belangrijk deze exacte oplossingen te weten wanneer de snelheid constant is. Er kan namelijk gekeken worden of de numerieke benadering bij constante snelheid wel deugt door deze te vergelijken met de exacte oplossing voor de golf. Voor beide numerieke en exacte methode kan er namelijk een simulatie gemaakt worden en dus het verloop van de golf bekeken worden. Als hier al de numerieke methode afwijkt van de exacte, zal er waarschijnlijk bij niet-constante snelheid c nog grotere fouten optreden bij het simuleren van de golf. Een belangrijke opmerking is dat deze oplossing alleen exact is voor constante c; als de snelheid variabel is kan hier geen gebruik van gemaakt worden. Er zal dus een numerieke methode geconstrueerd worden die accuraat de golf kan simuleren, ook bij variabele c.

3

F

OURIERTRANSFORMATIE

"Elke functie is te schrijven als een opsomming van sinussen en cosinussen", claimde Fourier in 1800. Hoewel hij toen niet geloofd werd, is het later toch uitgegroeid tot een zeer belangrijke ontdekking binnen de analyse. Zo is de Fourierreeks inderdaad een opsomming van sinussen en cosinussen die elke periodieke functie kan representeren. Door te kijken naar de frequenties van deze sinussen en cosinussen, kan je de periodieke functie weergeven in het frequentiedomein: als een sinus met bepaalde frequentie een grote amplitude heeft, zorgt dit voor een grote piek in het frequentiedomein op die bepaalde frequentie. We kunnen deze techniek uitbreiden tot de Fouriertransformatie, waarmee ook niet-periodieke functies bekeken kunnen worden op het frequentiedomein. Deze transformatie zal nog goed van toepassing komen bij het numeriek benaderen van de tweede afgeleide in de golfvergelijking. De theorie in dit hoofdstuk is voor-namelijk afkomstig van het boek van Olver[5]. Daarbij hielpen de notities van Osgood[10]waar de theorieën van Fourier en het frequentiedomein op intuïtieve wijze worden gepresenteerd.

3.1

Fourierreeksen

Een Fourierreeks is een reeks van cosinussen en sinussen die naar een periodieke extensie van een functie convergeert. Met periodieke extensie wordt bedoeld dat van een functie een interval genomen wordt en die achter elkaar geplakt wordt. Op deze manier verkrijgen we een periodieke functie die overeen komt met de originele functie op een zeker interval.

Definitie 3.1. Gegeven een stuksgewijs C1 functie f (x), dan kan er een 2l-periodieke extensie f (x) gemaakt worden. De fourierreeks die naar f (x) convergeert wordt gegeven door:

a0 2 + ∞ X j=1 ³ ajcos jπx l + bjsin jπx l ´ ,

waar de fouriercoëfficienten bepaald worden door:

a_j=1 l Z l −l f (x) cosjπx l dx, bj= 1 l Z l −l f (x) sin jπx l dx.

Een periodieke golf wordt dus opgesplitst in sinussen en cosinussen met verschillende fre-quenties. Deze verzameling frequenties noemen we het spectrum. Als een golf periode T heeft,

3.2. FOURIERTRANSFORMATIE

dan worden de frequenties gegeven door _Tn_{, n ∈ Z. Een golf is bandlimited als er een grootste} N is waar voor alle |n| > N, geldt dat cn= 0. Met de coëfficiënten cn die bij de verschillende

frequenties n_T horen kan je uiteindelijk een overzicht geven van de frequenties en hoe "zwaar" ze meewegen in de oorspronkelijke golf; door deze coëfficiënten wordt de periodieke functie in het frequentiedomein gerepresenteerd. We kunnen tot slot bovenstaande reeks omschrijven naar complexe e-machten, wat in het vervolg gebruikt zal worden. De reeks die dan ontstaat, wordt als volgt gegeven:[5]

(3.1) ∞ X j=−∞ cje −i jπx l , waar c_j= 1 2l Z l −l f (x)e−i jπxl dx.

3.2

Fouriertransformatie

Met de Fouriertransformatie kunnen we ook van aperiodieke functies het spectrum bepalen wat verderop goed van pas zal komen. Het idee achter deze transformatie is, is dat het limiet van de periode van de Fouriercoëfficiënten naar oneindig genomen wordt en zo een oneindige integraal krijgt. Als gevolg van de periodiciteit T die naar oneindig gaat, krijgen we nu ook een continuum aan frequenties, waar we voor de periodieke functie nog een discrete verzameling hadden. Dit spectrum laat wederom zien hoe zwaar elke frequentie meeweegt en is het resultaat als een Fouri-ertransformatie wordt toegepast. Met de afleiding achterwege, wordt de FouriFouri-ertransformatie gegeven door[5] ˆ f (k) =p1 2π Z _∞ −∞ f (x)e−ikx dx.

We kunnen deze functie zien als een lineaire transformatie, omdat het een functie gedefinieerd op het ruimtedomein afbeeldt op een functie gedefinieerd op het frequentiedomein en het voldoet aan de eisen van een lineaire transformatie. Een korte notatie voor deze operator wordt als volgt gegeven:

ˆ

f (k) = F [f (x)].

Op een zelfde manier kunnen we ook vanaf de functie op het frequentiedomein terug naar de functie op het ruimtedomein met de inverse Fouriertransformatie:

F−1_{[ ˆf (k)] = f (x) =} 1 p 2π Z ∞ −∞ ˆ f (k)eikxdk,

Deze inverse lineaire transformatie zorgt ervoor dat er van de functie op het frequentiedomein weer terug naar de functie op het ruimtelijke domein gegaan kan worden:

F−1£

F [f (x)]¤ = f (x).

3.2.1 Golfvergelijking transformeren

Met de transformatie kan de golf gerepresenteerd worden door zijn verschillende frequenties. Behalve dat dit vaak een stuk minder complex beeld van de golf geeft, kunnen we hier ook andere

3.2. FOURIERTRANSFORMATIE

handige dingen mee doen. Waar de transformatie in deze scriptie toegepast zal gaan worden, is de eendimensionale golfvergelijking:

∂2_u

∂t2(x, t) = c 2∂2u

∂x2(x, t).

We kunnen vergelijking omschrijven om deze in het frequentiedomein te bekijken, door te vermenigvuldigen met e−ikxen p1

2π en te integreren naar x, ofwel de Fouriertransformatie van x

toepassen. Laten we beginnen met de rechterhelft. Neem f (x, t) =∂u_∂x(x, t), na vermenigvuldiging en integratie zoals bovenstaand volgt:

1 p 2π Z _∞ −∞ ∂f ∂x(x, t)e −ikx _{dx =} 1 p 2π ³ f (x, t)e−ikx¯¯ ¯ ∞ −∞− Z _∞ −∞f (x, t)(−ik)e −ikx _dx´ =p1 2_πik Z ∞ −∞ f (x, t)e−ikx dx =p1 2_πik Z ∞ −∞ ∂u ∂x(x, t)e −ikx _dx =p1 2π ³

u(x, t)e−ikx¯¯ ¯ ∞ −∞− Z ∞ −∞u(x, t)(−ik)e −ikx _dx´ = (ik)2p1 2_π Z ∞ −∞

u(x, t)e−ikx dx = −k2u(k, t),ˆ

waar we aannemen dat u(±∞) =∂u_∂x(±∞) = 0. Voor de linkerkant is er een andere aanpak. Neem hiervoor dezelfde f waarop de Fouriertransformatie toegepast wordt, zodat het volgende afgeleid kan worden: 1 p 2_π Z ∞ −∞ ∂f ∂t(x, t)e −ikx _{dx =} 1 p 2_π Z ∞ −∞ lim h→0 f (x, t + h) − f (x, t) h e −ikx _dx =p1 2πh→0lim 1 h ³Z ∞ −∞f (x, t + h)e −ikx _{dx −}Z ∞ −∞ f (x, t)e−ikx dx´ =p1 2_πh→0lim p 2π_{f (k, t + h) −}ˆ p2πf (k, t)ˆ h =_∂∂ t ˆ f (k, t) = ∂ ∂t 1 p 2π Z ∞ −∞ f (x, t)e−ikx dx = ∂ ∂t 1 p 2π Z ∞ −∞ ∂u ∂t(x, t)e −ikx_dx = ∂ 2 ∂t2u(k, t),ˆ

waar in de laatste stap hetzelfde proces wordt afgegaan als met f . Op deze manier kunnen we de golfvergelijking bekijken in het frequentiedomein, wat resulteert in een gewone differentiaalver-gelijking:

(3.2) ∂

2

∂t2u(k, t) = −cˆ

3.2. FOURIERTRANSFORMATIE

3.2.2 Harmonische oplossing

Tot slot wordt een belangrijk resultaat gegeven die beschreven wordt in het boek van Cohen[4]. Hiervoor zullen we kijken naar de differentiaalvergelijking (3.2). De oplossing van deze differen-tiaalvergelijking is van de vorm

(3.3) _{u(k, t) = A(k)e}ˆ ickt_{+ B(k)e}−ickt.

Als op (3.2) ook de Fouriertransformatie van t wordt toegepast, kan de vergelijking omgeschreven worden tot:

−ω2ˆˆu(k,ω) = −c2k2ˆˆu(k,_ω).

Hieruit volgt de zogenaamde dispersierelatie ω= ck. Op deze manier is in te zien dat er een belangrijke relatie bestaat tussen de drie waardes, waar k aangeduid wordt als het golfgetal en_ω als het frequentiegetal. Nu is (3.3) te schrijven als

ˆ

u(k, t) = A(k)eiωt+ B(k)e−iωt.

Als vervolgens de inverse Fouriertransformatie van x hierop wordt toegepast, komt de volgende oplossing tot stand:

(3.4) _{u(x, t) =p}1 2π

Z ∞

−∞

A(k)ei(ωt+kx)_{+ B(k)e}i(−ωt+kx)dk.

Op deze wijze kan dus bij constante c elke golf geschreven worden als superpositie van har-monische golven ei(kx±ωt) met verschillende amplitudes. Eigenschappen van oplossingen van de golfvergelijking met constante snelheid kunnen dus bestudeerd worden door te slechts te kijken naar de oplossing ei(kx±ωt).

4

N

UMERIEKE METHODES

In de vorige hoofdstukken is de kennis opgedaan om een begin te kunnen maken aan het oplossen van het voorwaartse probleem. De bedoeling in dit hoofdstuk is dat er een numeriek schema gemaakt wordt om vanuit begincondities zo nauwkeurig mogelijk de golven te kunnen construeren die in het medium zullen propageren. Het numerieke schema zal verkregen worden door de tweede afgeleides van de golfvergelijking numeriek te benaderen met bepaalde technieken. Dit zal op een discreet domein in tijd en ruimte plaatsvinden waardoor nauwkeurigheid van het schema afhankelijk is van het aantal punten dat gebruikt wordt in deze domeinen; de punten in het tijdsdomein hebben een zekere onderlinge afstand, net zo voor de punten in het plaatsdomein. Op deze manier kan met behulp van de numerieke benadering van de golfvergelijking een formule gegeven worden die de drukwaardes op het volgende tijdsstip kan berekenen aan de hand van de drukwaardes op voorgaande tijdstippen. In de praktijk zal de voorkeur gaan naar waar relatief weinig punten gebruikt worden, opdat de rekentijd binnen de perken blijft.

Het numerieke schema k-space waar de toolboxK-WAVEgebruik van maakt, wordt in onder andere de artikelen van Treeby et al.[1]en Tabei et al.[2] beschreven en bestudeerd. Dit schema maakt deels gebruik van eindige differenties (finite differences). Eindige differenties is een bekende methode om afgeleides numeriek te benaderen. Er zijn verschillende varianten van eindige differenties: ze variëren in nauwkeurigheid bij de benadering van de afgeleides, maar ook in de rekentijd dat nodig is. De variant die bestudeerd zal worden, is de tweede orde centrale eindige differenties; de grootte van de fouten die deze variant oplevert is van tweede orde. Als eerst zal deze numerieke methode aan de hand van het boek van Cohen[4]en de webpagina van Grigoryan[18]worden bekeken en toegepast op de golfvergelijking. We zullen vervolgens aan de hand van de artikel van Tabei et al. een alternatief voor de benadering van de tweede afgeleide bekijken waar gebruik wordt gemaakt van Fouriertransformaties. Met deze twee technieken om tweede afgeleides te kunnen benaderen zal uiteindelijk de k-space methode volgen.

We zullen de numerieke methodes slechts in een dimensie bekijken; argumentatie door middel van berekeningen wordt anders te complex maar de methode komt wel in elke dimensie overeen. Het stelsel eerste orde vergelijkingen (2.2) en (2.3) en de tweede orde vergelijking (2.1) kunnen beide numeriek benaderd worden om een numeriek schema te krijgen voor het berekenen van de druk een volgend tijdstip. De methode die gepresenteerd wordt is gebaseerd op de tweede orde

4.1. TWEEDE ORDE EINDIGE DIFFERENTIES

vergelijking. Om een inzicht te geven hoe het eerste orde stelsel gebruikt kan worden om een numeriek schema te verkrijgen, heb ik deze met eerste orde eindige differenties uitgewerkt in de Appendix. De numerieke methodes op beide modellen kunnen worden uitgebreid naar twee en drie dimensies.[2]

4.1

Tweede orde eindige differenties

We zullen de tweede afgeleides van de eendimensionale golfvergelijking benaderen met eindige differenties. Voor het gemak beschrijven we nogmaals expliciet deze vergelijking:

∂2_u

∂t2(x, t) = c 2∂2u

∂x2(x, t).

Om te beginnen zullen we een raster definiëren waar de druk u op elk paar (xn, tm) gedefinieerd

zal worden. De rasterpunten worden voor de plaats gegeven door x_n+1= xn+∆x en evenzo voor

de tijd door t_{m+1= tm+}∆t. De rasterpunten (x_n, t_m) liggen dus in de tijdsrichting op afstand∆t en in de plaatsrichting op∆x van elkaar. Volgens de tweede orde centrale eindige differentie kan de tweede tijdsafgeleide op het punt (x, t) benaderd worden met de volgende vergelijking:[7]

∂2_u

∂t2(x, t) ≈

u(x, t +∆t)−2u(x, t)+ u(x, t −∆t) (∆t)2 .

Als dit ook voor de plaatsafgeleide gedaan wordt, volgt de golfvergelijking met numerieke benadering voor de tweede afgeleides:

u(x, t +∆t)−2u(x, t)+ u(x, t −∆t) (∆t)2 = c

2u(x +∆x, t)−2u(x, t)+ u(x −∆x, t)

(∆x)2 .

De schrijfwijze van de druk op de verschillende punten zal aangeduid worden met u(xn, tm) = umn.

Op deze manier zullen de tweede afgeleides op het rasterpunt (xn, tm) zoals bovenstaande

vergelijkingen worden aangeduid door:

∂2 ∂x2u m n = um n+1− 2u m n + um_n−1 (∆x)2 ∂2 ∂t2u m n = um+1_n − 2umn + um−1n (∆t)2

Met deze notatie volgt de numerieke benadering van de golvergelijking op het punt (xn, tm):

(4.1) u m+1 n − 2umn+ um−1n (∆t)2 = c 2u m n+1− 2umn + um_n−1 (∆x)2 .

Als we bovenstaande vergelijking omschrijven zien we dat de waardes {um_n}bepaald kan worden door het volgende tijdsstapschema te gebruiken:

(4.2) um+1_{n = s(u}m_{n+1+ u}m_{n−1) + (1 − s)2u}m_{n − u}m−1_n ,

4.1. TWEEDE ORDE EINDIGE DIFFERENTIES

x_n−1, t_m x_n, t_m x_n+1, t_m

x_n, t_m−1 xn, tm+1

Figuur 4.1: De druk op (xn, tm+1) wordt bepaald door drukwaardes op voorgaande tijdstippen

twee opeenvolgende rijen de waardes op de punten {(xn, tm−1)} en {(xn, tm)} voor alle n en vaste m

gegeven zijn, dan kan ook de druk op de rijen erboven benaderd worden. Hiermee kan dus steeds per tijdstip bepaald worden hoe groot de druk is op de verschillende posities en op deze manier wordt gesimuleerd hoe de golf zou propageren.

Merk echter wel op dat het hier om een benadering gaat! De methode van eindige differenties is niet exact en brengt dus een fout met zich mee. De grootte van de fout kan met Taylors ontwikkeling bepaald worden. We kunnen namelijk kijken naar Taylors ontwikkeling van de volgende twee functies:

u(x +∆x, t) = u(x, t)+∆x_∂∂ xu(x, t) + (∆x)2 2 ∂2 ∂x2u(x, t) + (∆x)3 6 ∂3 ∂x3u(x, t) + ..., u(x −∆x, t) = u(x, t)−∆x ∂ ∂xu(x, t) + (∆x)2 2 ∂2 ∂x2u(x, t) − (∆x)3 6 ∂3 ∂x3u(x, t) + ....

Tel deze twee bij elkaar op, dan volgt er na enig omschrijven: u(x +∆x, t)−2u(x, t)+ u(x −∆x, t)

(∆x)2 =

∂2

∂x2u(x, t) + O¡(∆x) 2_¢.

Ditzelfde geldt voor de tijdsafgeleide, waar dus de volgende vergelijking verkregen wordt: u(x, t +∆t)−2u(x, t)+ u(x, t −∆t)

(∆t)2 =

∂2

∂t2u(x, t) + O¡(∆t) 2_¢.

Deze benadering brengt bij beide afgeleides een fout met zich mee. Zo volgt er voor de numerieke benadering van de waarde um+1_n dat deze een fout van grootteO¡(∆x)2+ (∆t)2¢ heeft.

Hoewel dit een intuïtief goed te begrijpen idee is, is de uitvoering niet erg effectief wegens de relatief grote fout. Voor het simuleren van golven levert dit veel problemen op; de gesimuleerde golf kan sneller of langzamer propageren dan de werkelijke golf waardoor je fouten krijgt. Het is natuurlijk mogelijk om het raster dichter te maken, maar dat zal een stuk meer rekenkracht vereisen. Vandaar dat er andere methodes te vinden zijn om het proces van simuleren effectiever

4.2. k-SPACE

te maken. Zo kan er bijvoorbeeld gekozen worden voor een hogere orde schema of een aangepast schema om de accuraatheid te vergroten, waarvan enkele in het boek van Cohen gepresenteerd worden.[4] Voor deze scriptie hebben we echter genoeg aan dit schema aangezien deze niet complex is, maar toch goed te gebruiken is in het alternatieve schema dat gepresenteerd zal worden.

4.2

k-Space

We hebben net gekeken naar een methode die een numerieke benadering geeft voor de golf-vergelijking. Er zal nu een alternatief numeriek schema gepresenteerd worden dat k-space genoemd wordt omdat de plaatsafgeleides zal worden benaderd met behulp van Fouriertransfor-maties, waar het frequentiedomein (Fourierdomein/k-domein) een rol speelt. De tijdsafgeleide zal wederom met behulp van tweede orde eindige differenties benaderd worden. Ondanks de fout van de eindige differenties blijkt deze methode een stuk nauwkeuriger te zijn dan de vorige om twee redenen. Behalve dat de plaatsafgeleide met Fouriertransformaties exact opgelost kan worden, zal er bovendien nog een term geïntroduceerd worden in het schema dat de fout van de eindige differenties verkleint. We zullen zien dat deze methode zelfs exact is wanneer de geluidssnelheid in het medium constant is. Onder de aanname dat de eigenschappen van het medium maximaal 5% afwijken van die van water, heeft deze methode een grote accuraatheid voor relatief weinig rasterpunten. Wanneer er naar media gekeken wordt waar de geluidssnelheid grote sprongen maakt, zal eerder voor hogere orde eindige differenties bij de tijdsafgeleide gekozen worden.[1]

Ook bij het stelsel van eerste orde vergelijkingen kan k-space toegepast worden; de plaatsaf-geleides kunnen exact bepaald worden met Fouriertransformaties en dezelfde soort correctieterm kan in het schema geïntroduceerd.[2]Het voordeel van numerieke schema’s gebaseerd op het stelsel eerste orde vergelijkingen is dat er gebruik gemaakt kan worden van staggered grids (zie Appendix), wat de nauwkeurigheid van het schema nog meer vergroot. Tevens is voor het eerste orde stelsel gemakkelijk een Perfectly Matched Layer (PML) te implementeren, die er voor zorgt dat de golven op de randen van het computationele domein geabsorbeerd worden (en er dus geen reflecties ontstaan), wat in de praktijk van belang zal zijn. Bij de tweede orde vergelijking k-space methode is een dergelijke PML niet gemakkelijk in het schema te implementeren. Een reden om de k-space methode voor de tweede orde vergelijking te bestuderen is dat het concept minder complex en computationeel efficiënter is; het akoestische veld wordt namelijk alleen gedefiniëerd door de akoestische druk, in plaats van de akoestische druk en de snelheid van de deeltjes.[2]

4.2.1 Afleiding k-space methode

In dit deel zal bekeken worden hoe de benaderingen uitgevoerd worden en waar de correctieterm vandaan komt. Op deze manier kunnen we het tijdsstapschema van k-space afleiden. We zullen de afleiding maken vanuit de golfvergelijking met een constante c0. De golfvergelijking is dan

4.2. k-SPACE van de vorm: (4.3) ∂ 2_u ∂t2(x, t) = c 2 0 ∂2_u ∂x2(x, t).

Deze vergelijking zullen we weer numeriek benaderen op een raster met dezelfde onderlinge afstanden∆t en ∆x. Om te beginnen zal de tijdsafgeleide met eindige differentie te benaderd worden. Zoals we net zagen kan de tweede afgeleide benaderd worden door

∂2_u

∂t2(x, t) ≈

u(x, t +∆t)−2u(x, t)+ u(x, t −∆t) (∆t)2 .

De plaatsafgeleide zal op een totaal andere manier benaderd worden. Zoals in het hoofdstuk Fouriertransformaties getoond is, geldt F £∂2u

∂x2(x, t)¤ = −k2u(k, t) = −kˆ 2F [u(x, t)]. Hieruit volgt

dat dat het rechterlid van (4.3) te schrijven is als

∂2_u

∂x2(x, t) = −F

−1_£k2_{F [u(x, t)]¤.}

De numerieke benadering op het raster voor deze Fouriertransformaties kan op zeer accuraat uitgevoerd worden door gebruik te maken van de (Inverse) Fast Fourier Transform.[1] Het is een techniek die in 1965 door Cooley en Tukey gepubliceerd is.[21]De fout van de FFT bij een zekere afstand ∆x is haast verwaarloosbaar vergeleken met bijvoorbeeld de fout van tweede eindige differenties waar dezelfde afstand gebruikt zou worden. Omdat voorwaartse en inverse transformatie ingebouwde functies inMATLABzijn, zal ik me hier niet verder aan wijden. Op deze manier komt al een verbeterde numerieke benadering voor de golfvergelijking tot stand:

u(x, t +∆t)−2u(x, t)+ u(x, t −∆t) = −c20F−1£k2F [u(x, t)]¤.

Nu is er alleen nog het probleem dat de benadering voor de tijdsafgeleide een relatief grote fout met zich meebrengt. We zagen net namelijk al dat de benadering met eindige differenties een fout ter grootte vanO¡(∆t)2¢ met zich meebrengt. Zoals al eerder is vermeld, is het mogelijk naar eigenschappen te kijken van oplossingen door slechts te kijken naar de oplossing van de vorm ei(kx+ωt). Wat gebeurt er met deze oplossing als we hierop de tweede orde centrale eindige differentie toepassen?

ei(kx+ω(t+∆t))_{− 2e}i(kx+ωt)_{+ e}i(kx+ω(t−∆t)) (∆t)2 = e i(kx+ωt)eiω∆t− 2 + e−iω∆t (∆t)2 = ei(kx+ωt)2 cos(ω∆t)−2 (∆t)2 = −ω2ei(kx+ωt)sinc2(ω∆t 2 ).

Hier wordt bij de laatste gelijkheid eerst gebruik gemaakt van 2 − 2cos(ω∆t) = 4sin2(ω∆t₂ ) en vervolgens de term sinc(x) =sin(x)_x geintroduceerd. We zien dat de afgeleide niet exact is, maar er een extra term bij de afgeleide komt: sinc2(ω∆t₂ ). Als we deze wegdelen, krijgen we een dus

4.2. k-SPACE

gecorrigeerde eindige differenties waar de fout opgeheven is en dus een exacte waarde van de tweede afgeleide geeft.

Als we bovenstaande toepassen op de golfvergelijking, kan er een exacte benadering van de golfvergelijking gegeven worden in de volgende termen:

u(x, t +∆t)−2u(x, t)+ u(x, t −∆t) (∆t)2_sinc2₍ω∆t

2 )

= −c20F−1£k2F [u(x, t)]¤.

Met behulp van de dispersierelatie_{ω= ck kan bovenstaande vergelijking omgeschreven worden.} Hieruit volgt de volgende vergelijking, waar de gehele operatie die op u wordt toegepast aan de rechterkant de tweede orde k-space operator genoemd wordt:

(4.4) _{u(x, t +}∆t)−2u(x, t)+ u(x, t −∆t) = −c2F−1£k2_(∆t)2_sinc2₍crefk∆t

2 )F [u(x, t)]¤.

Voor de afleiding werd een constante c0gebruikt, maar die is nu vervangen door een c die variabel

kan zijn. Wat waarschijnlijk ook opvalt is dat er nu cref in de correctieterm staat. Hier moet

namelijk een juiste referentiesnelheid gekozen worden, die zo veel mogelijk in de buurt van c ligt, zodat de correctieterm daadwerkelijk de fout verkleint. Voor c0 constant wordt uiteraard

cref= c0 gekozen. Omdat crefconstant is kan deze niet gelijk aan c genomen worden als deze

laatste variabel is; dan moet er een passende waarde voor c_refgekozen worden.

Bovenstaande vorm kan tot slot omgeschreven worden naar een tijdsstapschema, waarmee het drukveld op een volgend tijdstip bepaald kan worden aan de hand van het drukveld op voorgaande tijdstippen. Voor het (x, t)-raster zoals eerder gedefiniëerd kan namelijk de verzameling van drukwaardes um+1_{= {u}m₁, u₂m, . . . , um_n}op het tijdstip t_m+1bepaald worden door:

(4.5) um+1_{= 2u}m_{− u}m−1_{− c}2₀F−1£k2_(∆t)2_sinc2₍crefk∆t

2 )F [u

m_]¤.

Let hier wel op het verschil van de drukwaardes tussen dit schema en het schema met eindige differenties (4.2). Bij eindige differenties kan de drukwaarde op het punt (xn, tm+1) bepaald

worden met behulp van waardes op enkele omliggende punten. Bij k-space wordt echter niet per punt, maar gelijk de gehele verzameling drukwaardes op het volgende tijdstip berekend vanuit de twee verzamelingen van drukwaardes op voorgaande tijdstippen.

t_m−1 tm

t_m+1

x1 x2 x3 . . . xn−1 xn

4.3. BEGINCONDITIES IMPLEMENTEREN

4.2.2 Vervolg k-space methode

Als we nu een oplossing invullen voor de golfvergelijking met constante snelheid, klopt ver-gelijking (4.4) dan? Om onszelf te overtuigen hiervan kunnen we dit direct nagaan door naar de oplossing u(x, t) = F(x − ct) + G(x + ct) te kijken, wat volgens (2.5) een algemene oplossing is. Om onderscheid te maken tussen Fouriertransformaties van tijd en plaats schrijven we respectievelijkF_t en F_x. Tevens kunnen we voor het gemak de volgende notatie gebruiken: Fx[u(x, t)] = Fxu(k, t). Met deze schrijfwijze komt er bij het rechterlid het volgende:

F [u(x, t)] =p1 2_π

Z ∞

−∞

u(x, t)e−ikx dx =p1 2π Z ∞ −∞F(x − c0 t)e−ikx _{dx +p}1 2π Z ∞ −∞G(x + c0 t)e−ikx dx =p1 2_π Z ∞ −∞ F(z)e−ik(z+c0t)_{d z +} 1 p 2_π Z ∞ −∞ G(z)e−ik(z−c0t) _{d z} = e−ikc0t 1 p 2π Z ∞ −∞

F(z)e−ikz _{d z + e}ikc0t 1

p 2π Z ∞ −∞ G(z)e−ikz d z = e−ikc0tF xF(k) + eikc0tFxG(k).

Na multiplicatie met onder andere de correctieterm moet er vervolgens weer terug getrans-formeerd worden. Als de equivalentie gebruikt wordt dat

(4.6) k2(∆t)2sinc2(c0k∆t 2 ) = 4 c2₀sin 2₍c0k∆t 2 ) = 1 c2₀(2 − 2cos(c0k∆t)) = 1 c2₀(2 − e ic0k∆t − e−ic0k∆t_),

kan met cref= c0 de inverse transformatie als volgt uitgeschreven worden:

F−1h_k2_(∆t)2_sinc2₍c0k∆t 2 ) ³ e−ikc0tF xF(k) + eikc0tFxG(k) ´i =−1 c2₀ ³

u(x, t +∆t)−2u(x, t)+ u(x, t −∆t)´.

De uitgebreide uitwerking van deze laatste gelijkheid staat in A2. Duidelijk is nu dat het rechterlid en linkerlid van (4.4) overeenkomt; de term met c0 valt weg tegen zichzelf in het

rechterlid die we nog niet meegenomen hebben in de berekening.

4.3

Begincondities implementeren

Er zijn twee methodes geïntroduceerd waarmee de drukwaardes op een volgend tijdstip benaderd kunnen worden wanneer de drukwaardes op de voorgaande twee tijdstippen bekend zijn. Hoe zit het dan op het begintijdstip t0? Op dit tijdstip is er geen drukwaardes bekend op het voorgaande

tijdstip, maar er zijn wel begincondities gegeven. Voor de eerste stap van t₀ naar t₁zullen de formules aangepast worden zodat de formules uitgedrukt worden in de begincondities. Stel dat

4.3. BEGINCONDITIES IMPLEMENTEREN

de beginwaardes als volgt gegeven zijn:

u(x, 0) = f (x),

∂

∂tu(x, 0) = g(x).

Voor de eerste stap op t0is de waarde u0ngegeven door f te discretiseren. Dat wil zeggen dat we

f bekijken op alle rasterpunt {x_n}, zodat volgt:

fn= f (xn) = u0n.

Net zo zal g gediscretiseerd worden door alleen naar de waardes op de rasterpunten {xn}te kijken.

Met behulp van eindige differenties voor eerste afgeleides kunnen we de waarde van g op elk rasterpunt (xn) schrijven als

g_{n= g(xn) =} ∂

∂tu(xn, 0) ≈

u1_{n− u}−1_n 2∆t .

Op deze manier kunnen de fictieve waardes u−1_n gesubstitueerd worden met behulp van de volgende vergelijking, die wegens de benadering wel een kleine fout met zich meebrengt:

u−1_{n = u}1_n− 2∆tg_n.

Voor de eerste tijdsstap kunnen deze waardes dus in de schema’s geïmplementeerd worden. Het schema van eindige differentie zal dan voor elk punt (xn, t1) de drukwaarde kunnen berekenen

met de formule:

(4.7) u1_n=s

2( fn−1+ fn+1) + (1 − s)fn+∆tgn.

Het schema van k-space zal met de volgende formule de hele verzameling drukwaardes berekenen op het tijdstip t1: (4.8) u1_{= f +}∆tg −1 2c 2¡ F−1_£k2_(∆t)2_sinc2¡ crefk∆t 2 ¢ F [f ]¤¢.

5

A

NALYSE VAN DE METHODES

Er zijn twee numerieke tijdsstapschema’s geïntroduceerd waarmee golven gesimuleerd kunnen worden. Er is echter al opgemerkt dat eindige differenties in een homogeen medium al een fout geven. Er kan een Gaussiaanse beginconditie u = e−x2 en _∂t∂_{u = 0 bekeken worden om het} verschil tussen eindige differenties en k-space globaal te kunnen aanschouwen. Deze kan met beide methodes voor een bepaald aantal tijdstappen gesimuleerd worden en vergeleken worden met d’Alemberts oplossing. Voor bijvoorbeeld∆x = 0,5 kunnen we uit deze figuren opmaken dat de methode van eindige differenties al zeer grote fouten oplevert, terwijl de k-space methode gelijk loopt met de exacte oplossing. Eindige differenties zal uiteraard nauwkeuriger worden

Eindige differenties k-space methode d’Alembert oplossing

Figuur 5.1: Voorbeeldsimulatie van een Gaussiaanse puls in het (u, x)-vlak na een bepaalde tijd op een raster met∆x = 0,5.

naarmate∆x kleiner wordt, maar daardoor wordt de rekentijd ook langer. Bovendien zien we al uit Figuur 5.1 dat k-space duidelijk nauwkeuriger is op hetzelfde raster. Om deze reden zal eindige differenties aan de kant geschoven worden en k-space beter bestudeerd worden.

Bij k-space moet er een referentiesnelheid gekozen worden. Wanneer deze gelijk aan de constante geluidssnelheid in het medium gekozen wordt, zal deze methode gelijk lopen met de exacte golf (voor een klein genoeg raster, zie volgend hoofdstuk). Maar wat gebeurt er wanneer de referentiesnelheid ongelijk aan de geluidssnelheid in het medium is? Als voorbeeld zijn de exacte golf en de k-space gesimuleerde golf na een bepaalde tijd in één figuur gezet voor∆x = 0,5 en is de golf die naar links propageert uitvergroot. In Figuur 5.2a zijn de twee golven niet van

5.1. HARMONISCHE GOLVEN

elkaar te onderscheiden. Hier is dan ook de referentiesnelheid gelijk aan de geluidssnelheid van het medium gekozen. In Figuur 5.2b is de golf wederom met k-space gesimuleerd, maar nu is de referentiesnelheid cref= 1000 m/s, terwijl de geluidssnelheid in het medium 1500 m/s is.

(a) Uitvergroting k-space met cref= c0 en

exacte oplossing. De golven lopen precies gelijk.

(b) Uitvergroting k-space met cref6= c0 en

exacte oplossing. Nu geeft k-space duidelijk een fout.

Figuur 5.2: Simulatie met k-space en de exacte oplossing in dezelfde figuren op een raster met ∆x = 0,5. Het is dus van belang om de referentiesnelheid goed te kiezen om (grote) fouten te voorkomen.

In de praktijk zal het medium waarin de golf zich voorplant niet homogeen zijn, omdat er naar biologisch weefsel gekeken zal worden. Er zullen dan delen zijn in het medium waar de geluidssnelheid niet overeenkomt met de gekozen referentiesnelheid. De correctieterm zal de fout dan niet geheel opheffen. We kunnen in kaart brengen hoeveel k-space zal verschillen van de werkelijke golf, door te kijken naar de fasefout (phase error). Bij dit type fout kijk je naar het faseverschil van de twee golven op een bepaald tijdstip. Om te bekijken hoe de nauwkeurigheid in termen van de fase van k-space verandert, kan er gekeken worden naar de exacte oplossing van een golf die met een vaste snelheid van 1500 m/s propageert en cref laten variëren. Een

handige manier om het verschil in fase te bekijken is door de golven op het frequentiedomein te bekijken. We zullen ons verdiepen in een resultaat uit het artikel van Treeby et al.[1]over de fout die zal ontstaan bij k-space wanneer crefongelijk aan c0is. In het volgende hoofdstuk zal de

rasterafstand∆x aan bod komen.

5.1

Harmonische golven

We zullen de analyse maken met behulp van de oplossing u(x, t) = ei(kx−ωt)+ ei(kx+ωt) die met een snelheid van c0= 1500 m/s propageren. Er kan dan gekeken worden naar het verschil in

5.1. HARMONISCHE GOLVEN

fase tussen de numerieke simulatie en de exacte oplossing na bijvoorbeeld 50 golflengtes. Als we golflengteλnemen, dan kunnen er enkele andere eigenschappen van de golf berekend worden:

1. De lengte van het interval waar exact 50 golflengtes in zitten: L = 50 ·λ. 2. De periode: T =1500λ .

3. Het frequentiegetal: kl=2_Lπ· 50.

4. Het golfgetal:ω_{l= 1500 · kl}

Voor de k-space methode is er tevens een geschikt raster nodig. Een verhouding tussen de ruimte-en tijdsafstand wordt gegevruimte-en door het zogruimte-enaamdeCFL-getal. Dit getal bepaalt op de volgende wijze de verhouding:

(5.1) CFL_{= c}2₀∆t

∆x

Met de waarde CFL_{= 0, 3 wordt een geschikte verhouding tussen de afstanden verkregen; dit}

geeft een goede balans tussen nauwkeurigheid en rekensnelheid in een zwak heterogeen medium. Verder zal er gebruik gemaakt worden van 4 rasterpunten per golflengte in de analyse van de fasefout.[1]

We kunnen de representatie van de harmonische golf op het frequentiedomein bekijken door de Fouriercoëfficiënten van de Fourierreeks te berekenen van u. Het deel van de oplossing met de harmonische golf die naar links beweegt wordt in het Fourierdomein gerepresenteerd wordt door:

cj(t) = 1 L Z L/2 −L/2 ei(klx+ωlt)_e−i j2πx/L _dx =1 Le iωlt Z L/2 −L/2 ei(kl− j2π/L)x_dx.

De laatste integraal in de vergelijking is wegens orthogonaliteit gelijk aan 0 als k_l6= j2_Lπ en gelijk aan 1 als k_l= j2_Lπ; ditzelfde geldt voor de andere harmonische golf. Oftewel, de enige coëfficient

ˆ

ujdie geen niet-nul oplossing is, is voor j = 50 en wordt gegeven door

ˆ uj(t) = 1 Le iωlt +1 Le −iωlt_.

We zullen voor het analyseren van het faseverschil tussen de originele oplossing en de numerieke oplossing kijken naar één richting van de golf, dus ei(kx+ωlt)_{. Deze wordt in het}

Fourierdomein gerepresenteerd door 1_Leiωlt_{. De golf beweegt naar links en daarbij zal duidelijk}

ook de oplossing in het Fourierdomein veranderd worden (want t verandert). Echter, na een bepaalde tijd zal de golf er weer hetzelfde uitzien omdat deze periodiek is; om precies te zijn na één periode T =2_ωπ. De oplossing ziet er dan in het Fourierdomein ook weer hetzelfde uit, want

eiωl(t+T)

= ei(ωlt+2π)