Dzień dobry mam drobny problem z metodą różnic skończonych. Nie potrawie zdefiniować warunku początkowego opisującego źródło położone w centrum rozważanego obszaru.
Pełna treść zadania prezentuje się następująco: Dane jest źródło gazu o stałej w czasie koncentracji (np. 80%). Współczynnik dyfuzji gazu wydostającego się ze źródła nie zależy od jego stężenia. Powietrze wokół źródła porusza się ze stałą prędkością (jednostajny, poziomy wiatr). Należy podać ustalony rozkład koncentracji gazu w powietrzu dla różnych stałych dyfuzji D, różnych prędkości wiatru v i różnych wydajności źródła. Należy zbadać dokładność rozwiązań. Rozważamy oczywiście wyłącznie "widok z góry" dzięki czemu problem jest dwuwymiarowy.
Równanie które muszę rozwiązać przedstawia się następująco:
nabla^2 u(x,y) = 0
Poniżej zamieszczam kod rozwiązujący owo równanie przy ustalonych warunkach brzegowych. Problem polega jedynie na tym że ustalając identyczny warunek na każdym brzegu nie otrzymuje informacji na temat źródła, które powinno być umieszczone w środku. // Rownanie: d2u/dx2+d2u/dy2=0 stacksize(100000000) clear(); clc
// Dlugosc boku L=1;
// Liczba elementów na jednym boku: nx=20;
// Temperatura na brzegu t=300;
// Warunki brzegowe: function[ub]=wbl(L,nx,indj) dx=L/nx; y=dx*(indj-1); ub=t; endfunction
function[ub]=wbp(L,nx,indj) dx=L/nx; y=dx*(indj-1); ub=t; endfunction
function[ub]=wbd(L,nx,indi) dx=L/nx; x=dx*(indi-1); ub=t; endfunction
function[ub]=wbg(L,nx,indi) dx=L/nx; x=dx*(indi-1); ub=t; endfunction // ---------------
// dlugosc elementu siatki dx=L/nx;
// Wspolczynniki rownania: a1=1/(dx^2); a2=-4/(dx^2);
// indeksy "i" i "j" for i=1:nx-1 for j=1:nx-1 indi((i-1)*(nx-1)+j)=j+1; indj((i-1)*(nx-1)+j)=i+1; end end
// Macierz wspolczynnikow A=zeros((nx-1)^2,(nx-1)^2); for i=1:(nx-1)^2 A(i,i)=a2; if indi(i)~=2 A(i,i-1)=a1; end
if indi(i)~=nx A(i,i+1)=a1; end
if indj(i)~=2 A(i,i-nx+1)=a1; end
if indj(i)~=nx A(i,i+nx-1)=a1; end end
// Wektor wyrazow wolnych B=zeros((nx-1)^2,1); for i=1:(nx-1)^2 if indi(i)==2 B(i,1)=B(i,1)-a1*wbl(L,nx,indj(i)); end
if indi(i)==nx B(i,1)=B(i,1)-a1*wbp(L,nx,indj(i)); end
if indj(i)==2 B(i,1)=B(i,1)-a1*wbd(L,nx,indi(i)); end
if indj(i)==nx B(i,1)=B(i,1)-a1*wbg(L,nx,indi(i)); end end
// Rozwiazanie ukladu rownan u=A^{-1}*B;
// Rozpisujemy wektor "u" w macierz for i=1:nx-1 for j=1:nx-1 upom(i,j)=u((i-1)*(nx-1)+j); end end
// Cala macierz U=zeros(nx+1,nx+1); for i=1:nx+1 U(1,i)=wbd(L,nx,i); U(nx+1,i)=wbg(L,nx,i); U(i,1)=wbl(L,nx,i); U(i,nx+1)=wbp(L,nx,i); end U(2:nx,2:nx)=upom;
// Wizualizacja x=[0:dx:1]; y=x;
xinit(); xset("colormap",hotcolormap(32)) plot3d1(x,y,U)
Z góry dziękuję za wszelką pomoc, ponieważ nie potrafię sobie poradzić z postawionym problemem od dłuższego czasu.