Praktická cvičení z Numerické matematiky

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

---

• Rovnice vedeni tepla, explicitní schema (Cely skript heateq.m)

du/dt = p d2u/dx2 + f(x,t)
Oblast: [a, b] x [ 0, T]

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

p = 0.25; 

fxt = 'exp( -t)';
u0  = 'sin(pi * x)';
alfa = '0 + 0 * t';
beta = '0 + 0 * t.^2';

Oblast Q_T:

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

Sit a volba kroku:

n = 10;
maxk = 100; 
h = (b - a) / n;
tau = T / maxk;

Stabilita:

sigma = p * tau / (h * h);

printf('Stabilita %g <= 0.5 ?\n', sigma);
if (sigma > 0.5)
    printf('\nNESPLNENO!\n\nVypocet 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 = eval(fxt);

Pocatecni a okrajove podminky:

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

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

Vlastni vypocet (cyklem):

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

Graf reseni:

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

---

• Rovnice vedeni tepla, explicitni schema (Program v jazyce C heateq.c)



---

• Vypocet nove casove vrstvy:

for(i = 1; i < n; i++)
{
  xi = a + i * h;
  ukplus1[i] = (1 - 2 * sigma) * uk[i] + sigma * uk[i - 1] 
                     + sigma * uk[i + 1] + tau * f(xi, tk);
}

• Zobrazeni reseni (ze souboru reseni.dat) v gnuplot:

splot "reseni.dat" matrix with lines linecolor 3 linewidth 0.1

% nebo zkracene

sp "reseni.dat" matrix w l lc 3 lw 0.1

---

• Rovnice vedeni tepla, explicitni schema, maticovy zapis (Cely skript zde)

Zmenime cast vlastni vypocet, nejprve pridame:
Pridame matici prechodu z jedne casove vrstvy na dalsi:

% Matice prechodu z jedne casove vrstvy na dalsi:

A = diag((1 - 2 * sigma) * ones(n - 1, 1));
A = A + diag(sigma * ones(n - 2, 1), -1);
A = A + diag(sigma * ones(n - 2, 1), 1);

% Vliv okrajových podmínek zahrneme do matice:

A = [ zeros(n - 1, 1), A, zeros(n - 1, 1) ];
A(1,1)      = -sigma;
A(end,end)  = -sigma;

Vlastni vypocet pak realizujeme takto:

for k = 1:maxk,
  newU = A * U(k,:)' + tau * f(k,2:end - 1)';
  U(k+1,2:end - 1) = newU';
end

---

• Rovnice vedeni tepla, implicitni schema, maticovy zapis (Cely skript zde)

Oproti maticovemu zapisu zmenime jen matici prechodu a vypocet Matice prechodu z jedne casove vrstvy na dalsi (implicitne, reseni soustavy rovnic):

A = diag((1 + 2 * sigma) * ones(n + 1, 1));
A = A - diag(sigma * ones(n, 1), -1);
A = A - diag(sigma * ones(n, 1), 1);

A(1,1:2)      = [ 1, 0];
A(end,end-1:end)  = [0, 1];

Vlastni vypocet:

for k = 1:maxk,
  newU = A \ (U(k,:)'  + tau * f(k,:)');
  U(k+1,:) = newU';
end

---