Rozwiążemy numerycznie równanie Laplace'a:
\(\displaystyle{ u_{xx}+u_{yy}= 0 \ \ (1) }\)
w obszarze:
\(\displaystyle{ R = \{(x,y): a < x < b, \ \ c < y < d \}}\)
z uwzgędnieniem warunków Dirichleta:
\(\displaystyle{ u(x,c) = f_{1}(x), \ \ u(x,d) = f_{2}(x), \ \ a \leq x \leq b, }\)
\(\displaystyle{ u(a,y) = g_{1}(y), \ \ u(b,y) = g_{2}(x), \ \ c \leq y \leq d. }\)
Zagadnieniem Laplace'a modelujemy na przykład rozkład temperatury w prostokątnej jednorodnej cienkościennej płytce, którą pokrywamy siatką punktów węzłowych:
\(\displaystyle{ x_{i} = a +ih, \ \ y_{j} = c +jk, \ \ i = 0,1,...,n, \ \ j= 0,1,..., m. }\)
Pochodne cząstkowe aproksymujemy różnicami centralnymi drugiego rzędu:
\(\displaystyle{ \frac{u_{i+1,j} -2u_{i,j} +u_{i-1,j}}{h^2} + \frac{u_{i,j+1} -2u_{i,j} +u_{ij-1}}{k^2} = 0 }\)
\(\displaystyle{ 2\left(\frac{h^2}{k^2} +1\right) - (u_{i-1,j} + u_{i+1,j}) - \frac{h^2}{k^2}\left(u_{i,j-1} +u_{i,j+1} \right) =0, }\)
\(\displaystyle{ i = 1,2,...,n-1, \ \ j=1,2,..., m-1 }\)
i warunki brzegowe:
\(\displaystyle{ u(x_{i},y_{0}) = f_{1}(x_{i}), \ \ i=1,2,...n-1,}\)
\(\displaystyle{ u(x_{i},y_{m}) = f_{2}(x_{i}), \ \ i=1,2,...n-1,}\)
\(\displaystyle{ u(x_{0},y_{j}) = g_{1}(y_{j}), \ \ j=1,2,...m-1,}\)
\(\displaystyle{ u(x_{n},y_{j}) = g_{2}(y_{j}), \ \ j=1,2,...m-1.}\)
Przykład
Rozwiążemy w programie MATLAB równanie Laplace'a:
\(\displaystyle{ u_{xx} + u_{yy} = 0, }\)
przyjmując następujące warunki brzegowe:
\(\displaystyle{ u(x,0) =0, \ \ u_(x,\pi) = \sin(x), \ \ 0\leq x \leq \pi, }\)
\(\displaystyle{ u(0,y) =0, \ \ u(\pi, y) = 0, \ \ 0 \leq y \leq \pi.}\)
Program w MATLAB
Kod: Zaznacz cały
function laplace(f1,f2,g1,g2,a,n,itmax,tol)
% Rozwiązanie równania Laplace'a w kwadracie z warunkami: u(x,0)=f1(x)
% u(x,a)=f2(x) and u(0,y)=g1(y) i u(a,y)=g2(y) używając pięciopunktowej
% metody różnicowej.
h=a/n;
z=0:h:a;
h
disp('_____________________________________________')
fprintf(' u= x\\y ')
fprintf('%4.2f ',z)
fprintf('\n')
disp('______________________________________________')
for i=1:n+1
u(i,1)=feval(f1,(i-1)*h);
u(i,n+1)=feval(f2,(i-1)*h);
u(1,i)=feval(g1,(i-1)*h);
u(n+1,i)=feval(g2,(i-1)*h);
end
iter=0;
err=tol+1;
while (err>tol)&(iter<itmax)
err=0;
for i=2:n
for j=2:n
oldu=u(i,j);
u(i,j)=(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1))/4;
res=abs(u(i,j)-oldu);
if (err<res);
err=res;
end
end
end
iter=iter+1;
end
for i=1:n+1
fprintf(' %4.2f',z(i))
for j=1:n+1
fprintf('%10.4f ',u(i,j))
end
fprintf('\n')
end
iter
Kod: Zaznacz cały
function f=f1(x)
f=0;
function f=f2(x)
f=sin(x);
function g=g1(x)
g=0;
function g=g2(y)
g=0;
>> laplace('f1','f2','g1','g2',pi,9,120,10^(-6))
h =
0.3491
________________________________________________________________________________________________________________________
u= x\y 0.00 0.35 0.70 1.05 1.40 1.75 2.09 2.44 2.79 3.14
________________________________________________________________________________________________________________________
0.00 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0.35 0.0000 0.0108 0.0228 0.0377 0.0570 0.0833 0.1196 0.1703 0.2416 0.3420
0.70 0.0000 0.0202 0.0429 0.0708 0.1072 0.1566 0.2248 0.3201 0.4541 0.6428
1.05 0.0000 0.0273 0.0578 0.0954 0.1444 0.2109 0.3028 0.4313 0.6118 0.8660
1.40 0.0000 0.0310 0.0658 0.1085 0.1643 0.2399 0.3444 0.4905 0.6957 0.9848
1.75 0.0000 0.0310 0.0658 0.1085 0.1643 0.2399 0.3444 0.4905 0.6957 0.9848
2.09 0.0000 0.0273 0.0578 0.0954 0.1444 0.2109 0.3028 0.4313 0.6118 0.8660
2.44 0.0000 0.0202 0.0429 0.0708 0.1072 0.1566 0.2248 0.3201 0.4541 0.6428
2.79 0.0000 0.0108 0.0228 0.0377 0.0570 0.0833 0.1196 0.1703 0.2416 0.3420
3.14 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
iter =
88