Praktická cvičení z Numerické matematiky

Numerické řešení PDR - 12. týden

• Vlnová rovnice, explicitní schéma (Cely skript waveeq.m)

Data ulohy: (uzivame vektorove operace, zapis operace "s teckou")

c2 = 4;  % (c na druhou)

fxt  = '0 * x';
u0   = 'sin(pi * x)';
u1   = '0';
alfa = '0';
beta = '0';

Oblast Q_T:

b = 1; a = 0; 
T = 10; 

Sit a volba kroku:

n = 32;
maxk = 800; 
h = (b - a) / n;
tau = T / maxk;

Stabilita:

sigma2 = c2 * (tau * tau) / (h * h);

printf('Stabilita: %g <= 1 ? \n', sigma2);
if (sigma2 > 1)
    printf('\nPOZOR: Vypocet nestabilni metodou!\n');
else
    printf('Splneno.\n');
end

Sit a zdrojovy clen:

xx = a:h:b; tt = 0:tau:T;

[xsit, tsit] = meshgrid(xx, tt);

x = xsit; t = tsit; f = 0 * x + eval(fxt);

Pocatecni a okrajove podminky:

U = 0 * xsit; % << reseni ma stejny rozmer jako sit

x = xx; U(1,:) = eval(u0);
U(2,:) = U(1,:) + tau * eval(u1); 

t = tt; U(:,1) = eval(alfa); U(:,end) = eval(beta);

Vlastni vypocet (cyklem):

for k = 2:maxk,
 for i = 2:(n-1)
  U(k + 1, i) = (2 - 2 * sigma2) * U(k, i) + sigma2 * U(k, i + 1) 
             + sigma2 * U(k, i - 1) - U(k - 1, i) + tau * f(k, i);
 end
end

Graf reseni:

  mesh(xsit, tsit, U);
% nebo 
  surf(xsit, tsit, U);

---

• Vlnová rovnice, implicitní schéma (Cely skript waveeq_imp.m)