//Феномен Рунге
function y=my_func(x)
y = 1./(1 + x.^2)
endfunction

function y=newton_method(xCount, x, xInterpolCount, Xn, Yn)
//заполняем нулями и высчитываем диагональ
    for i=1:xInterpolCount
        for j=1:xInterpolCount
            if i==j then
                M(i,j) = Yn(i);
            else
                M(i,j) = 0;
            end;
        end;
    end;

//высчитываем разделённые разности
    for i=1:+1:xInterpolCount-1
        for j=1:xInterpolCount-i
            M(j,j+i)=(M(j+1,j+i) - M(j,j+i-1)) / (Xn(j+i)-Xn(j));
        end
    end;

    for i = 1:xCount
        y(i) = Yn(1);
        if (xInterpolCount > 1)
        then
        p=1;
            for u=2:+1:xInterpolCount
                p = p*(x(i) - Xn(u-1));
                y(i) = y(i) + (M(1,u)*p);
            end
        end
    end
endfunction

function p=PolyLagr(m,n,x,xt,ft)
    // m кол-во Х  n- кол-во узлов
    //x-наши задаваемые Х     xt-узлы
    // ft - знач в узлах
    for k=1:m
        sum=0
        for i=1:n
            w(i)=1
            for j=1:n
                if (j ~= i)
            w (i) =w(i) * (x (k) - xt (j)) / (xt (i) - xt (j));
                end
            end
        sum = sum + ft(i) * w(i);
    end
    p(k) =sum;
    end
endfunction

function RP=RealP(n,f,f1)
    // f-данная функция   f1- интерполяц. полином
    for i=1:n
    RP(i)= abs(f(i)-f1(i))
    end
endfunction

// M- максимум в 10 производной
function Pr=Proizv(x,N) //n - степень
    i=0
   for  j=1:2:N //четные
       d1(j)=%i*((-1)^(i+1))*(factorial (j))./(1+%i*x)^(j+1)
       d2(j)=%i*((-1)^(j+i+1))*(factorial (j))./(1-%i*x)^(j+1)
       i=i+1
   end

   c=0
   for j=2:2:N //нечетное
       d1(j)=((-1)^(j+c))*(factorial (j))./(1+%i*x)^(j+1)
       d2(j)=((-1)^(j+c+1))*(factorial (j))./(1-%i*x)^(j+1)

     c=c+1
   end
   P=0
   for i=1:N
   P(i)=abs(0.5*d1(i)+0.5*d2(i))
end
Pr=P(N)
endfunction
//disp(Proizv(2,3))
function M=Max(N,X)
    for i=1:N
        f(i)=Proizv(X(i),N)
    end
    M=f(1)
    if N > 1 then
        for j=2:N
            if real(M) <= real(f(j))
                if imag(M)<=imag(f(j))

            M=f(j)
            end
        end
    end
    end
endfunction

function Ap=ApproximateP(nx,xx,nu,xu,M,N_Proiz) //приближенная оценка
    // nx- кол-во иксов
    //xx - исксы
    // ny - кол-во узлов
    // xy- узлы

    M=Max(N_Proiz,xx)
    //disp(M)
    for j=1:nx
        p=1
        for i=1:nu
            p=abs(p*(xx(j)-xu(i)))
        end

     Ap(j)=M*p/factorial(nx+1)
        end
endfunction

function Ch=Chebyshev(a, b, n) //узлы Чебышёва
for i=1:n
Ch(i)=0.5*(b+a)+0.5*(b-a)*cos((2*i+1)*%pi/(2*(n+1)));
end

endfunction


N_Proiz=10
xCount = 100;
xInterpolCount = 15;//кол-во узлов
a = -5;
b = 5;
h1 = (b-a)/(xCount - 1); //шаг для иксов функции
h2 = (b-a)/(xInterpolCount - 1); //шаг для узлов функции

x1 = [a:h1:b]; //иксы функции
x2 = [a:h2:b]; //узлы интерполяции

Ch=Chebyshev(a,b,xInterpolCount) //узлы интерпол по Чеб.

f1 = my_func(x1) //функция в иксах
F_Ch=my_func(Ch)
F1 = my_func(x2) //функция в узловых точках

fN = newton_method(xCount, x1, xInterpolCount, x2, F1);//полином Ньютона
fL=PolyLagr(xCount,xInterpolCount,x1,x2,F1)
RealN=RealP(xCount,f1,fN)
RealL=RealP(xCount,f1,fL)
Ap=ApproximateP(xCount,x1,xInterpolCount,x2)

fN_Ch = newton_method(xCount, x1, xInterpolCount, Ch, F_Ch);
fL_Ch=PolyLagr(xCount,xInterpolCount,x1,Ch,F_Ch)
Real_N_Ch=RealP(xCount,f1,fN_Ch)
Real_L=RealP(xCount,f1,fL_Ch)
Ap_Ch=ApproximateP(xCount,x1,xInterpolCount,Ch)

subplot(3,4,1)
plot(x1,f1,"r",x1,fN,"b");
xtitle('кр.-исходная функция, Син.-Ньютон')

subplot(3,4,2)
plot(x1,f1,"r",x1,fL,"g");
xtitle('кр.-исходная функция, зел.-Лагранж')

subplot(3,4,3)
plot(x1,RealN,"g");
xtitle('реал. погр. пол.Ньютона')

subplot(3,4,4)
plot(x1,RealL,'r');
xtitle('реал. погр. пол.Лагранжа')

subplot(3,4,5)
plot(x1,Ap,"r");
xtitle('приближенная оценка погрешности')

//Для Чебыш.

subplot(3,4,6)
plot(x1,f1,"r",x1,fN_Ch,"b");
xtitle('кр.-исходная функция, Син.-Ньютон')

subplot(3,4,7)
plot(x1,f1,"r",x1,fL_Ch,"g");
xtitle('кр.-исходная функция, зел.-Лагранж')

subplot(3,4,8)
plot(x1,Real_N_Ch,"g");
xtitle('реал. погр. пол.Ньютона')

subplot(3,4,9)
plot(x1,Real_L,'r');
xtitle('реал. погр. пол.Лагранжа')

subplot(3,4,10)
plot(x1,Ap_Ch,"r");
xtitle('приближенная оценка погрешности')