Condição de contorno periódica para a equação do calor em] 0,1 [

13

Vamos considerar uma condição inicial suave e a equação do calor em uma dimensão: no intervalo aberto , e vamos assumir que queremos resolvê-lo numericamente com diferenças finitas.

tu=xxu
]0,1[

Eu sei que para o meu problema a ser bem colocada, eu preciso dotá-lo de condições de contorno em e . Eu sei que Dirichlet ou Neumann funcionam bem.x=0x=1

Se eu tiver no primeiro caso pontos interiores para , então eu tenho desconhecidos: para , porque é prescrito nos limites.Nxk=kN+1k=1,,NNuk=u(xk)k=1,,Nu

No segundo caso, eu realmente tenho incógnitas e sei como usar o Neumann BC (homogêneo) para discretizar o laplaciano na fronteira, por exemplo, com a adição de dois pontos fictícios e e as igualdades:N+2u0,,uN+1x1xN+2

u1u12h=0=uN+2uN2h

Minha pergunta é sobre BC periódico. Eu tenho a sensação de que eu poderia usar uma equação, ou seja, mas talvez duas, e então usaria

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

mas eu não tenho certeza. Também não sei quantas incógnitas eu deveria ter. É ?N+1

bela83
fonte
Você tem condições de contorno de Dirichlet ou Neumann? O número de células fantasmas depende da ordem da aproximação das condições de contorno de Neumann.
ilciavo
@ilciavo, a pergunta é sobre condições de contorno periódicas.
Bill Barth

Respostas:

8

A melhor maneira de fazer isso é (como você disse), basta usar a definição de condições de contorno periódicas e configurar suas equações corretamente desde o início, usando o fato de que . De fato, ainda mais fortemente, as condições de contorno periódicas identificam com . Por esse motivo, você deve ter apenas um desses pontos no domínio da solução. Um intervalo aberto não faz sentido ao usar condições de limite periódicas, pois não há limite .u(0)=u(1)x=0x=1

Esse fato significa que você não deve colocar um ponto em pois é o mesmo que . Discretizando com pontos, você usa o fato de que, por definição, o ponto à esquerda de é e o ponto à direita de é .x=1x=0N+1x0 xNxN x0

esquema de uma grade periódica

Seu PDE pode ser discretizado no espaço como

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

Isso pode ser escrito em forma de matriz como que

tx=1Δx2Ax
A=[21001121000012110012].

Obviamente, não há necessidade de criar ou armazenar essa matriz. As diferenças finitas devem ser calculadas em tempo real, tendo o cuidado de lidar com o primeiro e o último pontos, especialmente conforme necessário.

Como um exemplo simples, o script MATLAB a seguir resolve com condições de contorno periódicas no domínio . A solução fabricada é usada, o que significa . Usei a discretização do tempo de Euler para simplificar e calculei a solução com e sem formar a matriz. Os resultados são mostrados abaixo.

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

Lote da condição inicial

Gráfico de solução em t = 0,5

Gráfico de solução em t = 1,0

Gráfico de solução em t = 2,0

Doug Lipinski
fonte
1
Ótima e simples solução !! no caso de alguém precisa dele aqui as implementações em Python
ilciavo
Ótimo ! Eu não queria ser muito específico sobre o método, pois eu mesmo utilizava a semi-discretização, por exemplo. Mas para a equação, você confirma que não há necessidade de especificar uma condição no derivativo? Mais precisamente, tal condição levaria a um problema incorreto? x
bela83
@ bela83 Você está certo de que não há necessidade de especificar nada além da condição inicial. Fazer isso resultaria em um sistema super determinado. Tudo o que você precisa fazer é ter um pouco de cuidado perto dos pontos finais do intervalo, para ter certeza de que as coisas são distribuídas periodicamente. Existem muitas maneiras válidas de fazer isso.
Doug Lipinski
-1

De acordo com isso, você deve impor condições de contorno periódicas como:

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

Uma maneira de discretizar a Equação de Calor implicitamente usando Euler para trás é

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

Resolvendo o sistema de equações

[Eu-ΔtΔx2UMA][você1n+1você1n+1vocêNn+1]=[você1nvocê2nvocêNn]

UMA=[-210 00 00 0...0 01-210 00 0...0 00 01-210 0...0 00 00 01-21...0 00 00 00 01-2...0 00 00 00 0...0 01-2]

você0 0vocêN+1

você1-vocêN=0 0você2-você0 02Δx-vocêN+1-vocêN-12Δx=0 0

vocêx

[0 010 00 0...0 0-10 0-10 010 0...10 0-10 00 00 00 00 0Eu-ΔtΔx2UMA0 00 00 00 00 00 00 0][você0 0n+1você1n+1você2n+1vocêNn+1vocêN+1n+1]=[0 00 0você1nvocê2nvocêNn]

O que lhe dá equações N + 2 e N + 2 incógnitas.

Você também pode se livrar da primeira das equações e das células fantasmas e chegar a um sistema de N equações e N desconhecidas.

ilciavo
fonte
você-1vocêNvocêN]0 0,1[xk=kN+1xN=NN+1
Nvocê0 0vocêN-1você-1vocêNvocê1vocêNvocê0 0vocêN+1
você0 0=vocêN+1N+1
vocêxvocê(0 0,t)=você(1,t)vocêx(0 0,t)=vocêx(1,t)
1
u1 = uN adiciona uma restrição adicional (equação u1-uN = 0) ao seu sistema
ilciavo