Kondisi batas untuk persamaan adveksi didiskritisasi dengan metode beda hingga

14

Saya mencoba mencari beberapa sumber untuk membantu menjelaskan bagaimana memilih kondisi batas saat menggunakan metode beda hingga untuk menyelesaikan PDE.

Buku-buku dan catatan yang saat ini saya memiliki akses untuk semua mengatakan hal serupa:

Aturan umum yang mengatur stabilitas di hadapan batas terlalu rumit untuk teks pengantar; mereka membutuhkan mesin matematika yang canggih

(A. Iserles Kursus Pertama dalam Analisis Numerik dari Persamaan Diferensial)

Misalnya, ketika mencoba menerapkan metode 2 langkah leapfrog untuk persamaan advection:

kamusayan+1=kamusayan-1+μ(kamusaya+1n-kamusaya-1n)

menggunakan MATLAB

M = 100; N = 100;

mu = 0.5;

c = [mu 0 -mu];
f = @(x)(exp(-100*(x-0.5).^2));

u  = zeros (M, N);
x = 1/(M+1) * (1:M);

u(:,1) = f(x);
u(:,2) = f(x + mu/(M+1));

for i = 3:N
    hold off;
    u(:,i) = conv(u(:,i-1),c,'same') + u(:,i-2);
    plot(x, u(:,i));
    axis( [ 0 1 0 2] )
    drawnow;
end

Solusinya berperilaku baik sampai mencapai batas, ketika itu tiba-tiba mulai berperilaku buruk.

Di mana saya bisa belajar bagaimana menangani kondisi batas seperti ini?

Simon Morris
sumber

Jawaban:

12

Tanggapan Sloede sangat teliti dan benar. Saya hanya ingin menambahkan beberapa poin untuk membuatnya lebih mudah dipahami.

Pada dasarnya, setiap persamaan gelombang memiliki kecepatan dan arah gelombang yang melekat. Untuk persamaan gelombang satu dimensi: kecepatan gelombang adalah konstan suatu yang menentukan tidak hanya kecepatan di mana informasi yang merambat di domain tetapi juga arahnya. Jika a > 0 , informasi bergerak dari kiri ke kanan dan jika a < 0 sebaliknya.

kamut+Sebuahkamux=0
SebuahSebuah>0Sebuah<0

Untuk metode lompatan katak, ketika Anda mendiskritasikan persamaan yang Anda dapatkan: atau: u n i =u n - 2 i +μ(u n - 1 i + 1 -u n - 1 i - 1 ) dimanaμ=-aΔt/Δx. Dalam kasus Anda,μ>0

kamusayan-kamusayan-22Δt+Sebuahkamusaya+1n-1-kamusaya-1n-12Δx=0
kamusayan=kamusayan-2+μ(kamusaya+1n-1-kamusaya-1n-1)
μ=-SebuahΔt/Δxμ>0yang diterjemahkan menjadi gelombang ke kiri. Sekarang jika Anda memikirkannya, gelombang yang bergerak ke kiri, hanya akan memerlukan kondisi batas di batas kanan karena semua nilai di sebelah kiri diperbarui melalui tetangga kanan mereka. Sebenarnya menentukan nilai apa pun di batas kiri tidak konsisten dengan sifat masalah. Dalam metode tertentu, seperti melawan angin sederhana, ini ditangani secara otomatis karena skema ini hanya melibatkan tetangga yang benar dalam stensilnya. Dalam metode lain, seperti katak leap, Anda harus menentukan beberapa nilai "benar".

Ini biasanya dilakukan melalui ekstrapolasi dari domain dalam untuk menemukan nilai yang hilang. Untuk masalah multi-dimensi dan non-kanonik, ini melibatkan menemukan semua vektor eigen dari fluks Jacobian untuk menentukan bagian mana dari batas yang benar-benar membutuhkan kondisi batas dan bagian mana yang memerlukan ekstrapolasi.

GradGuy
sumber
Secara fisik, apa artinya menggunakan persamaan ini dengan kondisi batas di sisi kiri dan kanan?
Frank
5

Jawaban umum
Masalah Anda adalah Anda tidak menyetel (atau bahkan menentukan) kondisi batas sama sekali - masalah numerik Anda tidak jelas.

Secara umum, ada dua cara yang mungkin untuk menentukan kondisi batas:

  1. Tetapkan kondisi batas dengan menentukan dan u 101kamu0kamu101 eksternal, misalnya melalui solusi yang tepat.
  2. Ubah stensil numerik sehingga hanya akan menggunakan informasi interior di perbatasan.

Ke mana Anda pergi sangat tergantung pada fisika masalah Anda. Untuk masalah tipe persamaan gelombang, orang biasanya menentukan nilai eigen dari fluks Jacobian untuk memutuskan apakah kondisi batas eksternal diperlukan, atau apakah solusi interior akan digunakan (metode ini biasa disebut 'upwinding').



kamusaya-1nkamusaya+1nsayan+1saya=1kamu0nkamu100n+1kamu101n

kamu1nkamu100n

Anda dapat menemukan versi modifikasi dari kode sumber Anda di bawah ini:

M = 100; N = 100;

mu = 0.5;

c = [mu 0 -mu];
f = @(x)(exp(-100*(x-0.5).^2));

u  = zeros (M, N);
x = 1/(M+1) * (1:M);

u(:,1) = f(x);
u(:,2) = f(x + mu/(M+1));

for i = 3:N
    hold off;
    %u(:,i) = conv(u(:,i-1),c,'same') + u(:,i-2);

    % Apply the numerical stencil to all interior points
    for j = 2:M-1
        u(j,i) = u(j,i-2) + mu*(u(j+1,i-1) - u(j-1,i-1));
    end

    % Set the boundary values by interpolating linearly from the interior
    u(1,i) = 2*u(2,i) - u(3,i);
    u(M,i) = 2*u(M-1,i) - u(M-2,i);

    plot(x, u(:,i));
    axis( [ 0 1 0 2] )
    drawnow;
end
Michael Schlottke-Lakemper
sumber
Jawaban yang bagus, dan selamat datang di scicomp, Sloede. Satu pertanyaan, biasanya saya melihat "angin" didefinisikan sebagai penggunaan stensil satu sisi di mana informasi diambil dari hanya satu batas domain. Apakah Anda bermaksud mengatakan itu dalam respons Anda?
Aron Ahmadia
1
Ya memang. Maaf, jika jawaban saya tidak cukup jelas. Akan tetapi, secara umum, "berliku" hanya berarti Anda mempertimbangkan arah aliran informasi. Itu tidak harus berarti bahwa Anda membuang satu sisi solusi sepenuhnya, itu hanya berarti Anda memberikan preferensi pada bagian dari solusi yang terletak pada arah "melawan arah angin".
Michael Schlottke-Lakemper
Jika Anda membuat N = 1000dan menjalankan kode sedikit lebih lama, Anda mendapati bahwa kode itu tidak berjalan sesuai harapan.
Simon Morris
Alasan untuk ini adalah bahwa solusi "perbaikan cepat" saya tidak sehat secara fisik, dan di atas itu agak sensitif terhadap osilasi palsu dalam solusi. Jangan gunakan ini untuk perhitungan ilmiah yang sebenarnya!
Michael Schlottke-Lakemper
2

Jadi saya telah melihat ini secara lebih rinci, dan tampaknya ini (setidaknya dalam kasus dasar yang saya tangani) tergantung pada kecepatan grup metode ini.

Metode leapfrog (misalnya) adalah:

kamusayan+1=kamusayan-1+μ(kamusaya+1n-kamusaya-1n)

Mencoba solusi formulir kamukn=esaya(ζkΔx+ω(ζ)nΔt) kami menemukan:

e2sayaωΔt=1+μesayaωΔt(esayaζΔx-e-sayaζΔx)

dosa(ωΔt)=μdosa(ζΔx)

dωdζ=cos(ζΔx)1-μ2ssayan2(ζΔx)[-1,1]

Sekarang kita perlu mengetahui kecepatan kelompok kondisi batas:

Metode saya :kamu1n+2=kamu1n+μkamu2n+1

Kita dapat menghitung kecepatan grup batas sebagai berikut:

2sayadosa(ωΔt)=μesayaζΔx

jadi untuk menemukan beberapa kecepatan grup yang batasnya memungkinkan kita perlu menemukan:

ω=cζ

cos(ζΔx)=0,μdosa(ζΔx)=2dosa(ζcΔt)

ζ=π2Δx akan memberi μ=2dosa(cμπ2) dimana solusi untuk c[-1,1]akan ada. (Untuk sebagian besar pilihanμ setidaknya)


Solusi yang saya temukan dalam literatur adalah mengambil kamu0n+1=kamu1n karena ini memiliki nomor gelombang batas yang terletak di luar [-1,1].

Saya masih sedikit lebih banyak untuk membaca tentang ini sebelum saya memahaminya sepenuhnya. Saya pikir kata-kata kunci yang saya cari adalah teori GKS.

Sumber untuk semua catatan A Iserles Part III ini


Perhitungan yang lebih jelas tentang apa yang telah saya lakukan dapat ditemukan di sini: http://people.maths.ox.ac.uk/trefethen/publication/PDF/1983_7.pdf

Simon Morris
sumber
-2

Kawan, saya sangat baru di situs ini. Mungkin ini bukan tempat untuk bertanya, tapi tolong maafkan saya karena saya sangat baru di sini :) Saya memiliki masalah yang sangat mirip, satu-satunya perbedaan adalah fungsi awal yang, dalam kasus saya, adalah gelombang kosinus. Kode saya adalah ini: hapus semua; clc; tutup semua;

M = 1000; N = 2100;

mu = 0,5;

c = [mu 0 -mu]; f = @ (x) 1- cos (20 * pi * x-0,025). ^ 2; u = nol (M, N); x = 0: (1 / M): 0,05; u (1: panjang (x), 1) = f (x); u (1: panjang (x), 2) = f (x - mu / (M)); x = linspace (0,1, M);

untuk i = 3: N tahan;

% Apply the numerical stencil to all interior points
for j = 2:M-1
    u(j,i) = u(j,i-2) - mu*(u(j+1,i-1) - u(j-1,i-1));
end

% Set the boundary values by interpolating linearly from the interior
u(M,i) =  2*u(M-1,i-1) - u(M-2,i-1);

plot (x, u (:, i)); sumbu ([0 1,5 -0,5 2]) digambar; % jeda akhir

Sudah ada kode ini di sini, tetapi untuk beberapa alasan, mungkin terkait dengan gelombang kosinus, kode saya gagal: / bantuan apa pun akan dihargai :) terima kasih!

John
sumber
2
Selamat datang di SciComp.SE! Anda harus membuat ini pertanyaan baru. (Jawaban hanya ditujukan untuk, sebenarnya, jawaban yang sebenarnya.) Jika Anda menggunakan "tanyakan tautan pertanyaan Anda sendiri" di bagian bawah (itu berwarna kuning tua dengan cahaya kuning, memang agak sulit untuk melihat apakah Anda tidak tahu itu ada di sana) , itu akan secara otomatis menautkan pertanyaan ke yang ini. (Anda juga dapat memasukkan tautan ke pertanyaan ini di dalam milik Anda.)
Christian Clason