KonuhovTask2

%a = const
Nx = 400;
X_L = 10;
U_0 = @(x) x/X_L*abs(sin(2*x/3));
mu = @(t) (1-exp(-t/10)+sin(2*t)^2/(1+t/100))/2;
Hx = X_L/Nx;
B = 1;
a = 5;
Ht = B*Hx/a;
gm = a*Ht/Hx;
x = zeros(1, Nx+1);
for i=0:Nx
    x(i+1)=i*Hx;
end
t_entered = 2;
Uj = zeros(1, Nx+1);
for i=0:Nx
    Uj(i+1)=U_0(x(i+1));
end
%точное решение для а>0
Uex = zeros(1, Nx+1);
for i=0:Nx
    if(x(i+1)>=a*t_entered)
        Uex(i+1)= U_0(x(i+1)-a*t_entered);
    else
        Uex(i+1)= mu(t_entered-x(i+1)/a);
    end
end
%приближенное решение для a>0
t = 0;
Uj1 = zeros(1, Nx+1);
while(t<=t_entered)
    Uj1(1)=mu(t);
    for i=1:Nx
        Uj1(i+1) = (1-gm)*Uj(i+1)+gm*Uj(i);
    end
    for i=0:Nx
        Uj(i+1)=Uj1(i+1);
    end
    t = t + Ht;
end
nexttile
plot(x, Uex)
hold on
plot(x, Uj, 'o')

%точно решение для a<0
Uex_flip = zeros(1, Nx+1);
for i=0:Nx
    if(X_L-x(i+1)>=abs(a)*t_entered)
        Uex_flip(i+1)= U_0(X_L-x(i+1)-abs(a)*t_entered);
    else
        Uex_flip(i+1)= mu(t_entered-(X_L-x(i+1))/abs(a));
    end
end
%приб решение для а<0
Uj_flip = flip(Uj);

nexttile
plot(x, Uex_flip)
hold on
plot(x, Uj_flip, 'o')

% a not const
Nx = 400;
X_L = 10;
U_0 = @(x) x/X_L*abs(sin(2*x/3));
mu = @(t) (1-exp(-t/10)+sin(2*t)^2/(1+t/100))/2;
Hx = X_L/Nx;
B = 1;
x = zeros(1, Nx+1);
for i=0:Nx
    x(i+1)=i*Hx;
end
t_entered = 2;
Uj = zeros(1, Nx+1);
for i=0:Nx
    Uj(i+1)=U_0(x(i+1));
end
Uj1 = zeros(1,Nx+1);
t = 0;
while(t<=t_entered)
    amax = sin(t/4.5-x(1)/2);
    for i=2:Nx+1
        a = sin(t/4.5-x(i)/2);
        if(abs(a)>amax)
            amax = abs(a);
        end
    end
    if(amax~=0)
        Ht=Hx/amax;
    else
        Ht=Hx;
    end
    Uj1(1)=mu(t);
    for i=0:Nx
        a(i+1) = sin(t/4.5-x(i+1)/2);
        gm = a(i+1)*Ht/Hx;
        nn=1;
        if(i==0)
            nn=0;
        end
        if(i==Nx)
            nn=2;
        end
        switch(nn)
            case 0
                if(gm<0)
                    Uj1(i+1)=(1+gm)*Uj(i+1)-gm*Uj(i+2);
                end
            case 1
                if(gm>=0)
                    Uj1(i+1)=(1-gm)*Uj(i+1)+gm*Uj(i);
                else
                    Uj1(i+1)=(1+gm)*Uj(i+1)-gm*Uj(i+2);
                end
            case 2
                if(gm>0)
                    Uj1(i+1)=(1-gm)*Uj(i+1)+gm*Uj(i);
                else
                    Uj1(i+1)=Uj(i+1);
                end
        end
    end
    for l=0:Nx
            Uj(l+1)=Uj1(l+1);
    end
    t = t+Ht;
end
nexttile
plot(x, Uj)

nexttile
plot(x, a)