Kondisi batas periodik untuk persamaan panas dalam] 0,1 [

13

Mari kita perhatikan kondisi awal yang halus dan persamaan panas dalam satu dimensi: dalam interval terbuka , dan mari kita asumsikan bahwa kita ingin menyelesaikannya secara numerik dengan perbedaan hingga.

tu=xxu
]0,1[

Saya tahu bahwa agar masalah saya dapat diposisikan dengan baik, saya harus memberinya syarat batas pada dan . Saya tahu bahwa Dirichlet atau Neumann bekerja dengan baik.x=0x=1

Jika saya memiliki dalam kasus pertama titik interior untuk , maka saya memiliki tidak diketahui: untuk , karena ditentukan pada batas-batas.Nxk=kN+1k=1,,NNuk=u(xk)k=1,,Nu

Dalam kasus kedua, saya benar-benar memiliki tidak diketahui , dan saya tahu cara menggunakan Neumann BC (homogen) untuk melenyapkan laplacian di perbatasan, misalnya dengan tambahan dua poin fiktif dan dan persamaannya:N+2u0,,kamuN+1x-1xN+2

kamu1-kamu-12h=0=kamuN+2-kamuN2h

Pertanyaan saya adalah tentang BC berkala. Saya merasa bahwa saya dapat menggunakan satu persamaan, yaitu tapi mungkin dua, dan kemudian saya akan menggunakan

kamu(0)=kamu(1)
xkamu(0)=xkamu(1)

tapi saya tidak yakin. Saya juga tidak tahu berapa banyak yang tidak diketahui yang harus saya miliki. Apakah ini ?N+1

bela83
sumber
Apakah Anda memiliki syarat batas Dirichlet atau Neumann? Jumlah sel hantu tergantung pada urutan perkiraan untuk Anda kondisi batas Neumann.
ilciavo
@ Silciavo, pertanyaannya adalah tentang kondisi batas periodik.
Bill Barth

Jawaban:

8

Cara terbaik untuk melakukan ini adalah (seperti yang Anda katakan) hanya menggunakan definisi kondisi batas periodik dan mengatur persamaan Anda dengan benar dari awal menggunakan fakta bahwa . Bahkan, lebih kuat lagi, kondisi batas periodik mengidentifikasi x = 0 dengan x = 1 . Karena alasan ini, Anda hanya boleh memiliki salah satu dari poin ini di domain solusi Anda. Interval terbuka tidak masuk akal ketika menggunakan kondisi batas periodik karena tidak ada batas .u(0)=u(1)x=0x=1

Fakta ini berarti Anda tidak boleh menempatkan titik pada karena sama dengan x = 0 . Diskritisasi dengan N + 1 poin, Anda kemudian menggunakan fakta bahwa menurut definisi, titik di sebelah kiri x 0 adalah x N dan titik di sebelah kanan x N adalah x 0 .x=1x=0N+1x0 xNxN x0

skema kotak periodik

PDE Anda kemudian dapat didiskresi dalam ruang sebagai

t[x0x1xN]=1Δx2[xN2x0+x1x02x1+x2xN12xN+x0]

Ini dapat ditulis dalam bentuk matriks sebagai mana A=[ - 2 1 0 0 1 1 - 2 1 0 0

tx=1Δx2Ax
A=[21001121000012110012].

Tentu saja tidak perlu untuk benar-benar membuat atau menyimpan matriks ini. Perbedaan yang terbatas harus dihitung dengan cepat, berhati-hati untuk menangani poin pertama dan terakhir secara khusus sesuai kebutuhan.

tu=xxu+b(t,x)
x[1,1)uRef(t,x)=exp(t)cos(5πx)b(t,x)=(25π21)exp(t)cos(5πx)
clear

% Solve: u_t = u_xx + b
% with periodic boundary conditions

% analytical solution:
uRef = @(t,x) exp(-t)*cos(5*pi*x);
b = @(t,x) (25*pi^2-1)*exp(-t)*cos(5*pi*x);

% grid
N = 30;
x(:,1) = linspace(-1,1,N+1);

% leave off 1 point so initial condition is periodic
% (doesn't have a duplicate point)
x(end) = [];
uWithMatrix = uRef(0,x);
uNoMatrix = uRef(0,x);

dx = diff(x(1:2));
dt = dx.^2/2;

%Iteration matrix:
e = ones(N,1);
A = spdiags([e -2*e e], -1:1, N, N);
A(N,1) = 1;
A(1,N) = 1;
A = A/dx^2;

%indices (left, center, right) for second order centered difference
iLeft = [numel(x), 1:numel(x)-1]';
iCenter = (1:numel(x))';
iRight = [2:numel(x), 1]';

%plot
figure(1)
clf
hold on
h0=plot(x,uRef(0,x),'k--','linewidth',2);
h1=plot(x,uWithMatrix);
h2=plot(x,uNoMatrix,'o');
ylim([-1.2, 1.2])
legend('Analytical solution','Matrix solution','Matrix-free solution')
ht = title(sprintf('Time t = %0.2f',0));
xlabel('x')
ylabel('u')
drawnow

for t = 0:dt:1
    uWithMatrix = uWithMatrix + dt*( A*uWithMatrix + b(t,x) );
    uNoMatrix = uNoMatrix + dt*(  ( uNoMatrix(iLeft) ...
                                - 2*uNoMatrix(iCenter) ...
                                  + uNoMatrix(iRight) )/dx^2 ...
                                + b(t,x) );
    set(h0,'ydata',uRef(t,x))
    set(h1,'ydata',uWithMatrix)
    set(h2,'ydata',uNoMatrix)
    set(ht,'String',sprintf('Time t = %0.2f',t))
    drawnow
end

Plot kondisi awal

Plot solusi pada t = 0,5

Plot solusi pada t = 1.0

Plot solusi pada t = 2.0

Doug Lipinski
sumber
1
Solusi hebat dan sederhana !! dalam kasus seseorang membutuhkannya di sini implementasi di Python
ilciavo
x
@ bela83 Anda benar bahwa tidak perlu menentukan lebih dari kondisi awal. Melakukannya akan menghasilkan sistem yang terlalu ditentukan. Yang perlu Anda lakukan adalah sedikit berhati-hati di dekat titik akhir interval untuk memastikan Anda membungkus dengan benar secara berkala. Ada banyak cara yang valid untuk melakukan ini.
Doug Lipinski
-1

Menurut ini, Anda harus memberlakukan ketentuan batas periodik sebagai:

u(0,t)=u(1,t)ux(0,t)=ux(1,t)

Salah satu cara mendiskritisasi Persamaan Panas secara implisit menggunakan Euler mundur adalah

un+1unΔt=ui+1n+12uin+1+ui+1n+1Δx2

Memecahkan sistem persamaan

[IΔtΔx2A][u1n+1u1n+1uNn+1]=[u1nu2nuNn]

A=[210000121000012100001210000120000012]

u0uN+1

u1uN=0u2u02ΔxuN+1uN12Δx=0

ux

[0100010101010100000IΔtΔx2A0000000][u0n+1u1n+1u2n+1uNn+1uN+1n+1]=[00u1nu2nuNn]

Yang memberi Anda persamaan N + 2 dan N + 2 tidak diketahui.

Anda juga dapat menyingkirkan yang pertama ke persamaan dan sel hantu dan tiba di sistem persamaan N dan N yang tidak diketahui.

ilciavo
sumber
Saya tidak mengerti pernyataan itu: "tambahkan dua sel tambahan kamu-1 dan kamuN"karena kamuN sudah menjadi titik masuk ]0,1[. Ada dalam pikiran sayaxk=kN+1 yang seperti itu xN=NN+1.
bela83
Ini hanya masalah pengindeksan. Anda mulai denganN sel (atau titik) dari kamu0 untuk kamuN-1 dan Anda menambahkan dua sel kamu-1 dan kamuN. Jika Anda memiliki sel darikamu1 untuk kamuN maka Anda perlu menambahkan sel di kamu0 dan kamuN+1
ilciavo
Baiklah kalau begitu saya tidak mengerti dua persamaan lagi! Yang pertama harus (dengan notasi dari pertanyaan):kamu0=kamuN+1Baik ? Tapi saya tidak mengerti yang kedua: mengapa tidak memilih perkiraan perbedaan terpusat? Akhirnya, itu membuatnyaN+1tidak diketahui, jika nilai pertama dan terakhir sama. Silakan bandingkan dengan situasi dengan (homogen) Neumann BC dalam pertanyaan.
bela83
Saya sudah mengubah notasi. Itu tergantung pada urutan aproksimasi untukkamux. Yang pertama berasalkamu(0,t)=kamu(1,t) dan yang kedua dari kamux(0,t)=kamux(1,t)
ilciavo
1
u1 = uN menambahkan batasan tambahan (persamaan u1 − uN = 0) ke sistem Anda
ilciavo