HMM_MIKRO_teil.m

function [u,w,index_n]=HMM_MIKRO_teil(T,u,epsilon,tau,dx,dt,n_x,M,index)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%        Berechnung der Loesung des mirkoskopischen Modells           %%%
%%%                      fuer einen Zeitschritt                         %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% function 
%
% Input:   T              Inverse Matrix von A
%          u              stueckweise konstante Rekonstruktion, AB
%          epsilon,tau    (Regularisierungs-) Parameter
%          dx, n_x        mikroskopische Ortsschrittweite, Anzahl
%                         mikroskopischer Ortsschritte
%          x              Stuetzstellen
%          dt, nt         mikroskopische Zeitsschritte, Anzahl
%                         mikroskopischer Zeitschritte
%          ul, ur         linker bzw. rechter Zustand
% Output:  u              Loesung des mirkoskopischen Modells auf dem mikro-
%                         skopischen Gitter 



if nargin==9
    index = sort(index);
    U = [];
    Ie = [];
    Ia = 1;
    I = [];
    W=[];
    for i=1:length(index)-1;
        if index(i+1)-index(i)==0
            error('Index kommt doppelt vor');
        elseif index(i+1)-index(i)>1     % Zusaetzliche Abfrage, dass mind 3 aufeinanderfolgende Werte
            Ia = [Ia, i+1];
            Ie = [Ie, i];
        end
    end
    Ie = [Ie,length(index)];
        
    for i=1:length(Ie)
        ni = Ie(i)-Ia(i)+1;  % Anzahl aufeinanderfolgender Stuetzstellen
        if ni<=2
           error('Zu wenig aufeinanderfolgende Stuetzstellen.') 
        end
        % Diskretisierungsmatrix A in IR^{(ni-2)x(ni-2)}
        n=ni-2;
        a = 2+dx^2/(epsilon^2*tau);
        c = epsilon^2*tau/(dx^2);
        A = c*(diag(a*ones(1,n))+diag(-ones(1,n-1),1)+diag(-ones(1,n-1),-1));
        A(1,1)=A(1,1)-1;    % evtl. BESSERE Randbedg. einfuegen!!!!!
        A(n,n)=A(n,n)-1;
        
        b = zeros(ni,1);
        for k=2:ni-1  
            b(k) = epsilon*(BL_D(u(index(Ia(i))+k-1))*(u(index(Ia(i))+k)-u(index(Ia(i))+k-1))-BL_D(u(index(Ia(i))+k-2))*(u(index(Ia(i))+k-1)-u(index(Ia(i))+k-2)))/(dx^2)-(BL_f(u(index(Ia(i))+k-1),M)-BL_f(u(index(Ia(i))+k-2),M))/dx;
        end 
        b(1) = [];
        b(end) = [];
        %%% Loesen des LGS Aw=b
        w = A\b;
        
        %%% Bestimmung der Loesung u mit Hilfe von w=u_t und oben berechnetem w
        u_local = u(index(Ia(i))+1:index(Ie(i))-1)+dt*w';
        U = [U,u_local];
        I = [I,index(Ia(i)+1):index(Ie(i)-1)];
        W = [W,w'];
    end
    u = U;
    w = W;
    index_n = I;
    
    
    
else
    b = zeros(n_x,1);
    for i=2:n_x-1 
        b(i) = epsilon*(BL_D(u(i))*(u(i+1)-u(i))-BL_D(u(i-1))*(u(i)-u(i-1)))/(dx^2)-(BL_f(u(i),M)-BL_f(u(i-1),M))/dx;
    end
    b(1) = [];
    b(end) = [];
    
    %%% Loesen des LGS Aw=b
    w = T*b;
    % Hinzunahme der RB w_1=0, w_nx=0
    w = [0;w;0];  
    %%% Bestimmung der Loesung u mit Hilfe von w=u_t und oben berechnetem w
    u = u+dt*w';    
end