Jak debugować kod numeryczny, co może być źródłem tego błędu oscylacyjnego?

Cichą wiedzę można uzyskać z doświadczenia, zastanawiałem się tylko, czy ktoś wcześniej widział coś podobnego do tego. Wykres pokazuje początkowy warunek (zielony) dla równania rada-dyfuzja, następnie roztwór przy iteracji 200 (niebieski), a następnie ponownie przy iteracji 400 (czerwony).

Rozwiązanie równania dyfuzyjno-doradczego pojawia się po kilku iteracjach. Liczba Pécleta i warunek CFL jest spełniony, , więc równania powinny być stabilne. Przewiduję, że mam błąd w kodzie numerycznym.$\mu\approx0.07$$C\approx 0.0015$

Tło. Dyskretyzacja stanowi zasadniczą różnicę zarówno w odniesieniu do warunków doradztwa, jak i rozpowszechniania. Uważam, że jest to pierwsza kolejność porad i druga kolejność rozpowszechniania. Zaimplementowałem to, stosując podejście o skończonej objętości (po raz pierwszy), w którym wartości współczynników (prędkości i współczynnika dyfuzji) na powierzchniach komórek określa się poprzez interpolację liniową z średnich komórek. Stosuję warunek brzegowy Robina na lewej i prawej powierzchni i ustawiam strumień na granicach na zero.

Jak debugujesz swój kod numeryczny? Czy ktoś wcześniej coś takiego zrobił, gdzie byłoby dobre miejsce, aby zacząć szukać?

Aktualizacja

• Oto moje osobiste notatki w stylu „książki laboratoryjnej” na temat implementacji metody o skończonej objętości do równania porada-dyfuzja, http://danieljfarrell.github.io/FVM/
• Kod źródłowy Python jest dostępny tutaj, http://github.com/danieljfarrell/FVM.git

Aktualizacja

Rozwiązanie nie może być prostsze! Właśnie popełniłem błąd znakowy dla terminu dyfuzji. To dziwne, jestem pewien, że nie opublikowałem tego, nie znalazłbym błędu! Jeśli ktoś chce podzielić się wskazówkami na temat debugowania swojego kodu numerycznego, nadal jestem zainteresowany. Nie mam metody, jest trochę trafiona i chybiam, próbuję różnych rzeczy, aby zdobyć wskazówki, ale ten proces może zająć tygodnie (czasem).

Dowód, że to działa (należy pamiętać, że przy metodzie o skończonej objętości wszystko, co musisz zrobić, aby obliczyć obszar, to suma szerokości wysokości dla wszystkich komórek, jeśli użyjesz metody integracji, takiej jak numpy.trapz, wyniki zawierają wartości liczbowe błąd metody trapezowej). Co tu się dzieje? Istnieje stała prędkość i współczynniki dyfuzji, ale przy zamkniętych warunkach brzegowych. Dlatego na granicy widzimy równowagę między polem prędkości popychającym w prawo a dyfuzją popychaną w lewo.$\times$

Zebrałem trochę mojego doświadczenia w debugowaniu kodów numerycznych tutaj: deal.II FAQ: debugowanie . Nie wiem, czy to pomogłoby ci w tym konkretnym przypadku, ale może w innych.