sanciyangtiao.m

根据给定函数的结点的函数值

n=length(x);
p=length(u);
%形成三转角方程组
for i=1:n-1
    h(i)=x(i+1)-x(i);
end
a2(1)=1;
g(1)=3*(y(2)-y(1))/h(1)-p1*h(1)/2;
for k=2:n-1
    a1(k-1)=h(k)/(h(k)+h(k-1));
    a2(k)=h(k-1)/(h(k)+h(k-1));
    g(k)=3*a2(k)*(y(k+1)-y(k))/h(k)+3*a1(k-1)*(y(k)-y(k-1))/h(k-1);
end
a1(n-1)=1;
g(n)=3*(y(n)-y(n-1))/h(n-1)-pn*h(n-1)/2;

%用第二章第六节追赶法求转三角方程组的解
b1(1)=2;
m(1)=g(1)/2;
b2(1)=a2(1)/b1(1);
for i=2:n
    b1(i)=2-a1(i-1)*b2(i-1);
    if(i~=n)
        b2(i)=a2(i)/b1(i);
    end
    m(i)=(g(i)-a1(i-1)*m(i-1))/b1(i);
end
for i=n-1:-1:1
    m(i)=m(i)-b2(i)*m(i+1);
end
%根据计算f(u(i))的近似值s(u(i)),期中x(k)在区间[x(k),x(K+1)]
for j=1:p
    for i=1:n
        if((u(j)>=x(i))&(u(j)<x(i+1)))
            k=i;
            break;
        end
    end
    s(j)=0;
    s(j)=s(j)+(h(k)+2*(u(j)-x(k)))*(u(j)-x(k+1))^2*y(k)/(h(k))^3;
    s(j)=s(j)+(h(k)-2*(u(j)-x(k+1)))*(u(j)-x(k))^2*y(k+1)/(h(k))^3;
    s(j)=s(j)+(u(j)-x(k))*(u(j)-x(k+1))^2*m(k)/(h(k))^2;
    s(j)=s(j)+(u(j)-x(k+1))*(u(j)-x(k))^2*m(k+1)/(h(k))^2;
end