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 QT:
    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)