task31konuhov

function [y] = progonka(N,hx,x,K,q_a,q_b,f)
    A = zeros(1,N);
    C = zeros(1,N+1);
    B = zeros(1,N);
    C(1)=1;
    for i=1:N-1
        A(i)=K(x(i+1)-hx/2);
        B(i+1)=K(x(i+1)+hx/2);
        C(i+1)=(A(i)+B(i+1));
    end
    A(N)=0;
    C(N+1)=1;
    F = zeros(1,N+1);
    F(1)=q_a;
    for i=2:N
        F(i)=hx^2*f(x(i));
    end
    F(N+1)=q_b;
    y=zeros(1,N+1);
    v1 = F(1)/C(1);
    v2 = F(N+1)/C(N+1);
    k1=B(1)/C(1);
    k2=A(N)/C(N+1);
    alpha = zeros(1,N);
    beta = zeros(1,N);
    alpha(1)=k1;
%     beta(1)=k2;
    beta(1)=v1;
    for i=2:N
        alpha(i)=B(i)/(C(i)-A(i-1)*alpha(i-1));
        beta(i)=(A(i-1)*beta(i-1)+F(i))/(C(i)-A(i-1)*alpha(i-1));
    end
    y(N+1)=(v2+k2*beta(N))/(1-k2*alpha(N));
    m=N;
    while(m>=1)
        y(m)=alpha(m)*y(m+1)+beta(m);
        m=m-1;
    end
end


%%%%%%%%%%%%%%%%%%%
format long
a=0; b=pi; q_a=0; q_b=1;
c1=(exp(-2*pi)-exp(-pi))./(2*(1-exp(-pi))); c2=c1+1;
u=@(x) exp(-x).*(-c1+sin(x)./2-cos(x)./2)-exp(-2*x)./2+c2;
f=@(x) exp(-x)+sin(x);
K=@(x) exp(x);
N=100;
hx=(b-a)/N;
x=zeros(1,N+1);
x(1)=a;
for i=1:N
    x(i+1)=x(i)+hx;
end
y = progonka(N,hx,x,K,q_a,q_b,f);
% plot(x,u(x),'o')
% hold on
% plot(x,y,'r')
% xlabel('x')
% ylabel('U(x),y(x)')
% legend('Точн. решение','Числен. решение')
% title('Погрешность:',norm(u(x)-y,inf))

Amplit = 0.05;
w = 100.5;
f_dtb = @(x) f(x)+Amplit*sin(w*x);

f2 = figure;
plot(x,f(x))
hold on
plot(x,f_dtb(x),'g')


Ampl = 0.05:-0.01:-0.05;
w = 200.5;
psi = zeros(1,numel(Ampl));
for i=1:numel(Ampl)
   y_dtb = progonka(N,hx,x,K,q_a,q_b,@(x) f(x)+Ampl(i)*sin(w*x));
   psi(i) = norm(y-y_dtb,inf);
end
f1 = figure;
plot(Ampl,psi);
xlabel('A')
ylabel('ψ')

% plot(x,f(x))
% hold on
% plot(x,f_dtb(x),'g')

% hold on
% plot(x,f_dtb(x),'r')
% f2 = figure;
% loglog([10^(-1) 10^(-2) 10^(-3)],[3.89*10^(-3) 3.91*10^(-5) 3.91*10^(-7)])
% xlabel('h')
% ylabel('ѱ')