Okresowe warunki brzegowe dla równania cieplnego w] 0,1 [

13

Rozważmy sprawne stan początkowy i równanie ciepła w jednym wymiarze: w otwartym przedziale ] 0 , 1 [ i załóżmy, że chcemy go rozwiązać numerycznie z różnic skończonych.

tu=xxu
]0,1[

Wiem, że aby mój problem został dobrze postawiony, muszę nadać mu warunki brzegowe przy i x = 1 . Wiem, że Dirichlet lub Neumann działają dobrze.x=0x=1

Jeśli mam w pierwszym przypadku punktów wewnętrznych x k = kN dlak=1,,N, a następnie mamNnieznanych:uk=u(xk)dlak=1,,N, ponieważujest określony na granicach.xk=kN+1k=1,,NNuk=u(xk)k=1,,Nu

W drugim przypadku naprawdę mam niewiadomych u 0 , , u N + 1 i wiem, jak używać (jednorodnego) Neumanna BC, aby zdyskretować laplaciana na granicy, na przykład przy połączeniu dwóch fikcyjnych punktów x - 1 oraz x N + 2 i równości:N+2u0,,uN+1x1xN+2

u1u12h=0=uN+2uN2h

Moje pytanie dotyczy okresowego pne. Mam wrażenie, że mógłbym użyć jednego równania, mianowicie ale może dwa, a następnie użyłbym ∂ x u ( 0 ) = x u ( 1 )

u(0)=u(1)
xu(0)=xu(1)

ale nie jestem pewien. Nie wiem też, ile niewiadomych powinienem mieć. Czy to ?N+1

bela83
źródło
Czy masz warunki brzegowe Dirichleta lub Neumanna? Liczba komórek-duchów zależy od kolejności przybliżenia warunków brzegowych Neumanna.
ilciavo
@ilciavo, kwestia dotyczy okresowych warunków brzegowych.
Bill Barth

Odpowiedzi:

8

Najlepszym sposobem na to jest (jak powiedziałeś) po prostu użycie definicji okresowych warunków brzegowych i prawidłowe ustawienie równań od samego początku przy użyciu faktu, że . W rzeczywistości, nawet silniej, okresowe warunki brzegowe utożsamiają x = 0 z x = 1 . Z tego powodu powinieneś mieć tylko jeden z tych punktów w domenie rozwiązania. Przedział otwarty nie ma sensu, gdy stosuje się okresowe warunki brzegowe, ponieważ nie ma granicy .u(0)=u(1)x=0x=1

x=1x=0N+1x0 xNxN x0

schemat okresowej siatki

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

tx=1Δx2Ax
A=[21001121000012110012].

Oczywiście nie ma potrzeby tworzenia ani przechowywania tej matrycy. Skończone różnice należy obliczać w locie, zwracając uwagę na pierwsze i ostatnie punkty specjalnie w razie potrzeby.

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

Wykres stanu początkowego

Wykres rozwiązania przy t = 0,5

Wykres rozwiązania przy t = 1,0

Wykres rozwiązania przy t = 2,0

Doug Lipiński
źródło
1
Świetne i proste rozwiązanie !! na wypadek, gdyby ktoś go tutaj potrzebował, implementacje w Pythonie
ilciavo
x
@ bela83 Masz rację, że nie trzeba określać niczego więcej niż warunek początkowy. Spowodowałoby to nadmiernie określony system. Wszystko, co musisz zrobić, to zachować ostrożność w pobliżu punktów końcowych interwału, aby mieć pewność, że okresowo odpowiednio się owijasz. Istnieje wiele ważnych sposobów, aby to zrobić.
Doug Lipinski,
-1

Zgodnie z tym należy nałożyć okresowe warunki brzegowe, ponieważ:

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

Jednym ze sposobów dyskretyzacji równania ciepła w sposób dorozumiany przy użyciu wstecznego Eulera jest

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

Rozwiązywanie układu równań

[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]

Co daje N + 2 równań i N + 2 nieznanych.

Możesz także pozbyć się pierwszych równań i komórek duchów i dojść do układu N równań i N nieznanych.

ilciavo
źródło
u1uNuN]0,1[xk=kN+1xN=NN+1
Nu0uN1u1uNu1uNu0uN+1
u0=uN+1N+1
uxu(0,t)=u(1,t)ux(0,t)=ux(1,t)
1
u1 = uN dodaje dodatkowe ograniczenie (równanie u1-uN = 0) do twojego systemu
ilciavo