Zajmiemy się numerycznym rozwiązaniem jednowymiarowego zagadnienia hiperbolicznego.
Rozpatrujemy równanie drgającej struny:
\(\displaystyle{ \alpha^2 u_{xx} = u_{tt}, \ \ 0< x < L, \ \ t>0,}\)
z warunkami brzegowymi:
\(\displaystyle{ u(0, t) = 0, \ \ u(L,t) =0, \ \ t >0, }\)
i warunkami początkowymi:
\(\displaystyle{ u(x,0) = f(x), \ \ 0 \leq x \leq L,}\)
\(\displaystyle{ u_{t}(x,0) = g(x), \ \ 0 \leq x \leq L.}\)
Uwzględniamy obszar prostokąta: \(\displaystyle{ R = \{(x,t): 0 \leq x \leq L, \ \ 0 \leq t \leq T \}, }\) który pokrywamy siatką \(\displaystyle{ (h,k) }\) i pochodne cząstkowe aproksymujemy różnicami centralnymi drugiego rzędu:
\(\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) }\)
\(\displaystyle{ u_{tt}(x_{i},t_{j}) = \frac{u(x_{i+1},t_{j}) -2u(x_{i},t_{j})+ u(x_{i-1},t_{j})}{k^2} + O(k^2) }\)
Podstawiamy różnice do równania \(\displaystyle{ (1) }\) opuszczamy składniki reszt i otrzymujemy równanie różnicowe:
\(\displaystyle{ \alpha^2 \frac{u_{j+1}^{j} -2u_{i}^{j}+u_{j-1}^{i} }{h^2}= \frac{u_{i}^{j-1} -2u_{i}^{j} +u_{i}^{j+1}}{k^2} \ \ (2) }\)
Kładąc \(\displaystyle{ \lambda = \alpha \frac{k}{h} }\) i przekształcając \(\displaystyle{ (2) }\) otrzymujemy schemat różnicowy trójpoziomowy
\(\displaystyle{ u_{i}^{j+1}= 2(1- \lambda^2)u_{i}^{j} +\lambda^2 (u_{j+1}^{j}+ u_{i-1}^{j}) -u_{i}^{j-1} \ \ (2)}\)
Można udowodnić , że ta metoda bezpośrednia (explicite) określona wzorem \(\displaystyle{ (2) }\) jest stabilnie numeryczna dla \(\displaystyle{ 0< \lambda \leq1 }\) i jej rząd dokładności wynosi \(\displaystyle{ O(k^3 +h^2 k^2).}\)
Jako przykład roziążemy równanie hiperboliczne:
\(\displaystyle{ 16 u_{xx} = u_{tt}, \ \ 0 < x < 1, \ \ t>0.}\)
spełniające warunki:
\(\displaystyle{ u(0,t) = u(1,t)=0, \ \ t>0,}\)
\(\displaystyle{ u(x,0) = \sin(\pi\cdot x), \ \ 0\leq x \leq 1,}\)
\(\displaystyle{ u_{t}(x,0) =0, \ \ 0\leq x \leq 1.}\)
Program w MATLAB
Kod: Zaznacz cały
function hiperbolic(f,g,c1,c2,L,T,h,k,alpha)
% Rozwiązanie równania hiperbolicznego z warunkami: u(x,0)=f(x), ut(x,0)=g(x)i
% u(0,t)=c1, u(L,t)=c2 trzypoziomowym schematem explicite
n=L/h;
m=T/k;
lambda=alpha*k/h;
z=0:h:L;
disp('___________________________________________')
fprintf(' t x = ')
fprintf('%4.2f ',z)
fprintf('\n')
disp('___________________________________________')
fprintf('% 4.2f ',0)
for i=1:n+1
u0(i)=feval(f,(i-1)*h);
fprintf('%10.6f ',u0(i))
end
fprintf('\n')
fprintf('% 4.2f %10.6f ',k,c1)
for i=1:n-1
u1(i+1)=(1-lambda^2)*feval(f,i*h)+lambda^2/2*(feval(f,(i+1)*h)...
+feval(f,(i-1)*h))+k*feval(g,i*h);
fprintf('%10.6f ',u1(i+1))
end
fprintf('%10.6f ',c2)
fprintf('\n')
for j=2:m
t=j*k;
fprintf('% 4.2f ',t)
u1(1)=c1; ui(1)=c1; u1(n+1)=c2; ui(n+1)=c2;
fprintf('%10.6f ',ui(1))
for i=2:n
ui(i)=2*(1-lambda^2)*u1(i)+lambda^2*(u1(i+1)+u1(i-1))-u0(i);
fprintf('%10.6f ',ui(i))
end
fprintf('%10.6f ',ui(n+1))
fprintf('\n')
u0=u1;
u1=ui;
end
Kod: Zaznacz cały
function f=f1(x,t)
f=sin(pi*x);
function g=g1(x,t)
g=0;
>> hiperbolic('f1','g1',0,0,1,0.5,0.2,0.05,4)
_______________________________________________________________________
t\ x = 0.00 0.20 0.40 0.60 0.80 1.00
_______________________________________________________________________
0.00 0.000000 0.587785 0.951057 0.951057 0.587785 0.000000
0.05 0.000000 0.475528 0.769421 0.769421 0.475528 0.000000
0.10 0.000000 0.181636 0.293893 0.293893 0.181636 0.000000
0.15 0.000000 -0.181636 -0.293893 -0.293893 -0.181636 0.000000
0.20 0.000000 -0.475528 -0.769421 -0.769421 -0.475528 0.000000
0.25 0.000000 -0.587785 -0.951057 -0.951057 -0.587785 0.000000
0.30 0.000000 -0.475528 -0.769421 -0.769421 -0.475528 0.000000
0.35 0.000000 -0.181636 -0.293893 -0.293893 -0.181636 0.000000
0.40 0.000000 0.181636 0.293893 0.293893 0.181636 0.000000
0.45 0.000000 0.475528 0.769421 0.769421 0.475528 0.000000
0.50 0.000000 0.587785 0.951057 0.951057 0.587785 0.000000