Zajmiemy się analizą numeryczną równania różniczkowego cząstkowego:
\(\displaystyle{ \alpha u_{xx} = u_{t}, \ \ 0 < x < L, \ \ 0 < t < T, \ \ (1) }\)
z warunkami brzegowymi:
\(\displaystyle{ u(0,t) = 0, \ \ u(L,t) =0, \ \ 0 \leq t \leq T }\)
i warunkiem początkowym:
\(\displaystyle{ u(x,0) = f(x), \ \ 0 \leq x \leq L .}\)
Na przykład modelujemy zjawisko rozchodzenia się ciepła w jednorodnym pręcie o długości \(\displaystyle{ L, }\) funkcja \(\displaystyle{ u(x, t) }\) oznacza wartość temperatury w chwili \(\displaystyle{ t }\) na długości \(\displaystyle{ x }\) pręta.
Zastosujemy metodę bezpośrednią (explicite) polegającą na przybliżeniu pochodnych cząstkowych w równaniu \(\displaystyle{ (1) }\) odpowiednimi różnicami centralnymi.
W tym celu budujemy siatkę punktów: \(\displaystyle{ x_{i}= h\cdot i, \ \ t_{j}= k\cdot j , \ \ i =0,1,..., n, \ \ j=0,1,..., m.}\)
\(\displaystyle{ u_{xx}(x_{i},t_{j}) = \frac{u(x_{i+1}, t_{j}) - 2u(x_{i},t_{j})+u(x_{i-1},t_{j})}{h^2} + O(h^2) , \ \ (2) }\)
\(\displaystyle{ u_{t}(x_{i}, t_{j}) = \frac{u(x_{i},t_{j+1}) - u(x_{i},t_{j})}{k} + O(k) \ \ (3) }\)
Uwzględniamy równania \(\displaystyle{ (2), (3) }\) w równaniu \(\displaystyle{ (1) }\)
\(\displaystyle{ \frac{u_{i}^{j+1} - u_{i}^{j}}{k} = \alpha \frac{u_{i+1}^{j}-2u^{j}_{i}+u_{i-1}^{j}}{h^2} \ \ (4) }\)
Podstawiamy \(\displaystyle{ \lambda = \alpha \frac{k}{h^2} }\) i rozwiązujemy równanie \(\displaystyle{ (4) }\) względem \(\displaystyle{ u_{i}^{j+1}. }\)
Otrzymujemy schemat różnicowy:
\(\displaystyle{ u_{i}^{j+1} = (1 -2\lambda)u_{i}^{j} + \lambda( u_{i+1}^{j} +u_{i-1}^{j}), \ \ i=1,2,...,n-1, \ \ j = 0,1, \ \ ... m-1.}\)
Można pokazać, że metoda explicite ma rząd dokładności \(\displaystyle{ O(k + h^2).}\)
Dla przykładu rozwiążemy równanie
\(\displaystyle{ u_{xx} = u_{t}, \ \ 0 < x < 1, \ \ , t >0 , }\)
z warunkami brzegowymi:
\(\displaystyle{ u(x,0) = \sin(\pi x) , \ \ 0\leq x \leq 1, \ \ u(0,t) = u(1,t) =0, \ \ t\geq 0.}\)
Program w MATLAB
Kod: Zaznacz cały
function rownanie_ciepla(f,c1,c2,L,T,h,k,alpha)
% Rozwiązanie równania ciepła z warunkami: u(x,0)=f(x)i u(0,t)=c1
% u(L,t)=c2 metodą różnic centralnych.
n=L/h; m=T/k;
lambda=alpha^2*k/(h^2);
z=0:h:L;
disp('________________________________________________________________')
fprintf(' t x = ')
fprintf('%4.2f ',z)
fprintf('\n')
disp('________________________________________________________________')
fprintf('% 5.4f ',0)
for i=1:n+1
u(i)=feval(f,(i-1)*h);
fprintf('%10.6f ',u(i))
end
fprintf('\n')
for j=1:m
t=j*k;
fprintf('% 5.4f ',t)
for i=1:n+1
if (i==1)
y(i)=c1;
elseif (i==n+1)
y(i)=c2;
else
y(i)=(1-2*lambda)*u(i)+lambda*(u(i+1)+u(i-1));
end;
fprintf('%10.6f ',y(i))
end;
fprintf('\n')
u=y;
end
Kod: Zaznacz cały
function f=f1(x,t);
f= sin(pi*x);
________________________________________________________________________________________________________________________________
t \ x = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
________________________________________________________________________________________________________________________________
0.0000 0.000000 0.309017 0.587785 0.809017 0.951057 1.000000 0.951057 0.809017 0.587785 0.309017 0.000000
0.0025 0.000000 0.301455 0.573401 0.789219 0.927783 0.975528 0.927783 0.789219 0.573401 0.301455 0.000000
0.0050 0.000000 0.294078 0.559369 0.769905 0.905078 0.951655 0.905078 0.769905 0.559369 0.294078 0.000000
0.0075 0.000000 0.286881 0.545680 0.751064 0.882929 0.928367 0.882929 0.751064 0.545680 0.286881 0.000000
0.0100 0.000000 0.279861 0.532327 0.732685 0.861322 0.905648 0.861322 0.732685 0.532327 0.279861 0.000000
0.0125 0.000000 0.273012 0.519300 0.714755 0.840244 0.883485 0.840244 0.714755 0.519300 0.273012 0.000000
0.0150 0.000000 0.266331 0.506591 0.697263 0.819682 0.861865 0.819682 0.697263 0.506591 0.266331 0.000000
0.0175 0.000000 0.259813 0.494194 0.680200 0.799623 0.840773 0.799623 0.680200 0.494194 0.259813 0.000000
0.0200 0.000000 0.253455 0.482100 0.663554 0.780055 0.820198 0.780055 0.663554 0.482100 0.253455 0.000000
0.0225 0.000000 0.247253 0.470303 0.647316 0.760966 0.800127 0.760966 0.647316 0.470303 0.247253 0.000000
0.0250 0.000000 0.241202 0.458793 0.631475 0.742343 0.780546 0.742343 0.631475 0.458793 0.241202 0.000000