Konuhov2zadanie

%aconst.m
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;
% U_0 = @(x) sqrt(x/X_L)*sqrt(abs(sin(4*x/3)));
% mu = @(t) 1/3*(1-exp(-t/15))*(cos(t)+sin(2*t))^2;
U_0 = @(x) (1-x/X_L)^5*abs(sin(x/2));
mu = @(t) exp(-t/20)*sqrt(abs(sin(t/5)));
Hx = X_L/Nx;
B = 1;
a = 1;
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 = 5;
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
f1 = figure;
f1.Position= [0 0 800 350];
plot(x, Uex)
hold on
plot(x, Uj, 'o')
legend('точ. решение', 'приб. решение')
xlabel('x')
ylabel('U(x,t)')

%точно решение для 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
% f = figure;
% f.Position = [0 0 800 350];
% plot(x, Uex_flip)
% hold on
% plot(x, Uj_flip, 'o')
% legend('точ. решение', 'приб. решение')
% xlabel('x')
% ylabel('U(x,t)')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%anotconst.m
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;
U_0 = @(x) sqrt(x/X_L)*sqrt(abs(sin(4*x/3)));
mu = @(t) 1/3*(1-exp(-t/15))*(cos(t)+sin(2*t))^2;
Hx = X_L/Nx;
B = 1;
x = zeros(1, Nx+1);
for i=0:Nx
    x(i+1)=i*Hx;
end
t_entered = 5;
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=B*Hx/amax;
    else
        Ht=B*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
f1 = figure;
plot(x, Uj)
xlabel('x')
ylabel('U(x,t)')

f2 = figure;
plot(x, a, 'r')
xlabel('x')
ylabel('a(x,t)')

% f1 = figure;
% f1.Position = [0 0 800 350];
% plot(x, Uj)
% hold on

% B=1.4;
% t = 0;
% Uj = zeros(1, Nx+1);
% for i=0:Nx
%     Uj(i+1)=U_0(x(i+1));
% end
% Uj1 = zeros(1,Nx+1);
% 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=B*Hx/amax;
%     else
%         Ht=B*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
%
% plot(x, Uj, 'o')
% legend('b=1','b=1.4')