Multigrid, PDE, programmering - boundary condition?

Hej

Vi håller på med ett projekt i en programmeringskurs där ska lösa en partiell differentialekvation. Man skulle kunna sammanfatta programmet på det här sättet:

En matris initieras med ettor på alla kanter, alltså längs med alla sidor, och nollor överallt annars.

För varje iteration(relaxation) går man igenom varje punkt och punktens värde är ett medelvärde av de kringliggande punkterna, med tillräckligt många iterationer kommer det till slut vara ettor överallt.

Uppgifterna i projektet är att tillämpa olika algoritmer för att få detta att fungera, bland annat jakobi, red/black, SOR och multigrid.

I de sista uppgifterna skall man applicera multigridmetoden. Vi har spenderat ungefär 6-7 timmar på bara att hitta information om detta och försöka förstå oss på den, och sedan lika mycket till att försöka översätta det till kod.

Programmet fungerar såhär just nu:

1. initiera matrisen
2. kör några iterationer med valfri annan relaxtionmetod
3. förminska matrisen (make it more coarse), gör detta 4 gånger och gör punkt 2 mellan varje gång.
4. förstora matrisen (interpolate the matrix, making it more fine), gör detta 4 gånger och gör punkt 2 mellan varje gång.

Matrisen förminskas genom att varannan punkt(både i x och y-led) blir en ny punkt i en annan matris:
Punkten i den nya består alltså av hälften av punkten i den gamla och 1/8 av dennes grannar.

Ungefär på samma sätt förstoras matrisen igen:
I den större(more finer) matrisen så kommer varannan punkt i båda leden vara samma som i den gamla(more coarser) och dess grannar kommer tilldelas hälften eller en fjärdedel av värdet punkten(the coarser). Men eftersom det är varannan punkt så kommer de ha gemensamma grannar, grannarna får då hälften av ena gamla(more coarser) punkten av hälften av nästa gamla punkt. Hörnpunkterna är gemensam granne för fyra punkter(coarser) och de kommer få en fjärdedel av varje sådan punkt. Det här är så krångligt, vet inte hur jag ska förklara det bättre...

Problemet som uppstår är att istället för att det är ettor överallt går det mot noll på kanterna och med tillräckligt många iterationer blir även punkterna i mitten noll.

Vi är väldigt osäkra huruvida man ska göra med "boundary condition", med de tidigare algoritmerna så har vi ju inte ändrat på storleken av matrisen. Nu när vi gör det, ska vi då lägga på ettor runt varje matris när vi förminskar den? Eller hur fungerar egentligen med "Dirichlet boundary condition"?

Görs det här i matlab eller är det riktig programmering? När man faltar (vilket jag tolkar det som att ni gör) i matlab kan man i sista argumentet ange hur den ska tolka punkter utanför matrisen när en kernel bara delvis överlappar den. Så om det nu görs i matlab borde ni snabbt kunna gå igenom de olika alternativen och se vilket som fungerar (om något). Jag tycker det är väldigt otydligt hur matrisen förminskas, först antyder du att den bara samplas: "Matrisen förminskas genom att varannan punkt(både i x och y-led) blir en ny punkt i en annan matris" för att sedan antyda att det är någon typ av bilinjär interpolation: "Punkten i den nya består alltså av hälften av punkten i den gamla och 1/8 av dennes grannar." Hur ligger det till?

Men det är så att jag går en kurs med huvudsakligen matematiker (datateknik själv) där alla uppgifter ska göras i matlab. Så efter att ha suttit i en timme för att få rätt på min trippelnästade for-loop (var en matris som representerade punkter i rummet) med ett dussin if-satser så kollar man på matematikernas lösning som har gjorts med en enda kodrad där de använt någon matlab-metod jag aldrig hört talas om, som dessutom gör att exekveringen går 10 ggr snabbare.

Detta är "riktig programmering" i C.

Jag förstår att det är otydligt, fattar knappt själv vad jag skrev. Ska försöka förklara det lite bättre.

Det du säger om bilinjär interpolation stämmer säkert, men eftersom vi inte har tagit någon kurs om partiella differentialekvationer så kan vi inga termer.


X: större matris
Y: mindre matris

När X förminskas till Y väljs punkter ut från X, hälften av denna punkten går till punkten i Y plus 1/8 av X grannar(ej "hörngrannar").


Så om en matris skulle se ut så här:
och man väljer ut de nya punkterna markerade med [] så här blir första punkten i Y lika med 0/2+(1+0+0+1)*1/8=1/4.

Det vi undrar egentligen är om man ska välja punkterna så här istället:

Problemet då är att en punkt på kanten eller i ett hörn inte har några grannar i en eller två riktningar.

Den första verkar ju vettig, man man behöver inte oroa sig för att kärnan hamnar utanför matrisen och den reducerar en 4x4 till en 2x2 matris, men du hävdar att det inte ger en korrekt lösning?

Vårat program är rekursivt, så det kan hända att något händer som inte får hända. Vi måste nog kolla igenom det ordentligt, kanske är det inte fel på hur vi reducerar matrisen utan hur vi sedan förstorar upp den igen? För någonstans blir det fel iallafall Vi drog iväg ett mail till lärarna i kursen idag, men han som är huvudsaklig ansvarig för kursen och förmodligen den som kan mest är utomlands för tillfället. Hoppas vi kan få svar i början på nästa vecka.

Kursen vi läser heter "Programmering av parallella system", som namnet antyder kommer vi dessutom behöva göra programmet parallellt senare. Men det är annan femma det får vi ta när det sekventiella fungerar som det ska helt enkelt.

Vad är det för något ni simulerar? Är det en sorts värmeledningsekvation ni löser? Varför förväntar ni er att få ettor överallt i matrisen (förhoppningsvis en ledande fråga)? Om ni behåller ettor i kanterna hela tiden (alltså att erat randvilkor, eller boundary condition, är att funktionsvärdet ska vara = 1) så kommer ni att simulera värmeledning där temperaturen utanför alltid är = 1 och temperaturen innanför = 0 från början. Eftersom temperatur tenderar att jämna ut sig bör allting gå mot 1.

Jag vet däremot inte hur ni utför alla steg, det verkar ju som att ni inte håller erat randvilkor om ni säger att det går mot noll även på kanterna. Hur gör ni för att räkna ut det nya värdet i kanterna? Det ni vill ha är som sagt ett randvilkor som säger att funktionsvärdet i kanterna hela tiden = 1, alltså ska ni inte räkna ut vad funktionsvärdet ska bli där utan helt enkelt bara sätta det till 1.