proj3.m

function fn = intcompare(x, f, a, b)

    hp = build_hermite(x, hermitecalc(x, f));
    sp = build_spline( x,splinecalc(x,f) );

    i_f = integral(f,a,b);
    sp_i = spline_integral(x, a, b, sp);
    hp_i = integral( matlabFunction(hp),a,b );
    
    first_row = [i_f,sp_i,hp_i];
    second_row = [ 0, abs(i_f-sp_i), abs(i_f-hp_i) ];
    third_row = [ 0, abs((i_f-sp_i)/(i_f)), abs((i_f-hp_i)/(i_f)) ];
    
    fn = [first_row;second_row;third_row];
end

function s = splinecalc(x, f)
    n = length(x);
    a=[];b=[];c=[];d=[];h=[];k=[];
    for i = 1:n
        a(i) = f( x(i) );
    end
    for j = 1:(n-1)
        h(j) = x(j+1) - x(j);
    end
    l = [1];
    m=[0];
    z=[0];
    for j = 2:(n-1)
        k(j) = ( a(j+1)-a(j) )*3/h(j) - ( a(j)-a(j-1) )*3/h(j-1);
        l(j) = 2*( x(j+1)-x(j-1))-h(j-1)*m(j-1);
        m(j) = h(j)/l(j);
        z(j) = ( k(j)-h(j-1)*z(j-1) )/l(j);
    end
    l(n) = 1;
    z(n) = 0;
    c(n) = 0;
    for j=(n-1):-1:1
        c(j) = z(j)-m(j)*c(j+1);
        b(j) = (a(j+1)-a(j))/h(j) - h(j)*( c(j+1)+2*c(j) )/3;
        d(j) = (c(j+1)-c(j))/(3*h(j));
    end
    for j=1:n-1
        s(1,j) = a(j);
        s(2,j) = b(j);
        s(3,j) = c(j);
        s(4,j) = d(j);
    end
end

function h = hermitecalc(x, f)
    n = length(x)-1;
    q = [];
    z=[];
    coef = [];
    for i=1:n+1
        z(2*i) = x(i);
        z(2*i-1) = x(i);
        q(2*i,1) = f( x(i) );
        q(2*i-1,1) = f( x(i) );
        
        syms k;
        y = diff( f(k) );
        g = matlabFunction(y);
        q(2*i,2) = g( x(i) );
        
         if i~=1
             q(2*i-1,2) = ( q(2*i-1,1)-q(2*i-2,1) )/( z(2*i-1)-z(2*i-2) );
         end
    end
    
     for i=3:(2*n+2)
         for j=3:i
              q(i,j) = ( q(i,j-1)-q(i-1,j-1) )/( z(i)-z(i-j+1) );
         end
     end
     
     for i=1:2*n+2
        for j=1:2*n+2
            if i==j
                 coef = [coef, q(i,j)];   
            end
        end 
     end
     h=coef;
end

function hp = build_hermite(pts, a)
    double_pts = repelem(pts,2);
    n = length(double_pts);
    syms p(x);
    p(x) = 0;
    for i=2:n
        tmp=1;
        for j=i-1:-1:1
            tmp = tmp*(x-double_pts(j));
        end
        p = p + a(i)*tmp;
         
    end
    p = p + a(1);
    hp = vpa(p,5);
end

function sp = build_spline(x, coef)
    n = length(x)-1;
    syms q(w);
    p(1) = q;
    for i=1:n
        syms polyn(m);
        polyn(m) = 0;
        for j=0:3
           if j==0
               polyn(m) = coef(1,i);
           else
               polyn(m) = polyn(m) + coef(j+1,1)*((m-x(i))^j);
           end
        end
        p(i) = polyn;
    end
    sp = p;
end

function int = spline_integral(x, a, b, sp)
    n = length(x);
    xx = x;
    int = 0;
    for i=1:n-1
        if x(i)==a
            break;
        end
        if x(i+1) > a
            xx = [x(1:i) a x(i+1:end)];
            break;
        end
    end
    for i=length(xx):-1:1
        if xx(i)==b
            break;
        end
        if x(i-1) < b
            xx = [xx(1:i) b xx(i+1:end)];
            break;
        end
    end
    
    for i=1:n-1
        if xx(i+1) > a
            int = integral( matlabFunction(sp(i)),a,xx(i+1) );
            idx = i;
            for j=i+1:n-1
                if j==idx && xx(j)==b
                    break;
                end
                if xx(j+1)==b
                    int = int + integral( matlabFunction(sp(j-1)),xx(j),b );
                else
                    int = int + integral( matlabFunction(sp(j-1)),xx(j),xx(j+1) );
                end
            end
        end
    end
end