
% Smisena uloha pro vlnovou rovnici         
% d^2u/dt^2 = c^2 d^2u/dx^2+f(x,t) ... v (a,b)x(0,T)              
% pocatecni:                       
% poloha  : u(x,0) = phi(x)       v <a,b>              
% rychlost: du(x,0)/dt = psi(x)  v <a,b> 
% okrajove
% u(a,t) = alfa(t) v <0,T>              
% u(b,t) = beta(t) v <0,T>              
%                                                      
%  METODA SITI - IMPLICITNI SCHEMA            

% vstupy:
% -------
 c = 1.5;        % POZOR! v zadani rovnice je c^2!!! 
               % kladna konstanta, rychlost sireni vlny
 f = @(x,t)3*x - 2*t;   % funkce na prave strane vlnove rovnice 
 a = 0;   % interval promenne x         
 b = 4; 
 %
 phi = @(x)2*x + 2*x.^2;    % okamzita pocatecni vychylka 
 psi = @(x)2*x  ;        % okamzita pocatecni rychlost 
 %
 A0 = phi(a)
 B0 = phi(b)
 A1 = psi(a)
 B1 = psi(b)
 alfa = @(t) A0 + A1*t;   % okrajova podminka v bode a
 beta = @(t) B0 + B1*t;   % okrajova podminka v bode b

% podminky souhlasu,t.j. 
% phi(a)== alfa(0),phi(b)==beta(0),
% psi(a)== d(alfa(0))/dt, psi(b)==d(beta(0))/dt,

h   = 0.4;    % prostorovy krok, 
tau = 0.05; % casovy krok,

sigma = c * tau / h;
sigma2 = sigma^2;

 
if sigma <= 1 
   v=['Podminka stability pro zadane h =',num2str(h),...
       ', tau=',num2str(tau),...
       ' a  c=',num2str(c),' je splnena. sigma =',num2str(sigma),...
       ', sigma^2=',num2str(sigma2)];
    disp(v);
else
   v=['Podminka stability pro zadane h=',num2str(h),...
         ', tau=',num2str(tau),...
      ' a c =',num2str(c),' neni splnena! sigma = ',num2str(sigma)];
   disp(v);
end

taumax = h/c

% výpocet
% -------
% body xi
x = a : h :b+h/10 %ověřit, zda x(end) je opravdu b
pocet_bodu=length(x);
% nulta vrstva
U0 = phi(x)
% prvni vrstva
U1 = phi(x) + tau * psi(x)
U1(1)   = alfa(tau); 
U1(end) = beta(tau);
% druha vrstva
U2(1) = alfa(2*tau);
U2(pocet_bodu) = beta(2*tau);
tau2=tau*tau;
for i = 2:pocet_bodu-1
    U2(i)=sigma2*(U1(i-1)+U1(i+1))+2*(1-sigma2)*U1(i)-U0(i)+tau2*f(x(i),tau);
end

U2





  
