mathanalysis304.m

数值分析的实验程序

N=10;
i=0:N;
x=-1+2/N*i;
fx=1./(1+25*(x.^2));
xx=-1:0.001:1;
ffx=zeros(1,length(xx));
for i=1:N
    ffx(find(xx==x(i)))=fx(i);
end
h=(max(x)-min(x))/N;
H=4*h*diag(ones(N+1,1))+h*diag(ones(N,1),1)+h*diag(ones(N,1),-1);
H(1,1)=H(1,1)/2;
H(N+1,N+1)=H(N+1,N+1)/2;
for i=1:N
    fi(i)=(fx(i+1)-fx(i))/h;
end
D=zeros(N-1,1);
for i=1:N-1
    D(i)=fi(i+1)-fi(i);
end
ddfx=D/h;
D=[0 
    D
    0];
D=6*D;
M=mathanalysis202(H,D);
for i=1:N
    F=find(xx>x(i)&xx<x(i+1));
    temp1=M(i)/6*((x(i+1)-xx(F))).^3;
    temp2=M(i+1)/6*(xx(F)-x(i)).^3;
    temp3=(h*h/6)*(M(i)+(M(i+1)-M(i))*(xx(F)-x(i))/h);
    ffx(F)=(temp1+temp2)/h+fx(i)+fi(i)*(xx(F)-x(i))-temp3;           
end
figure;
plot(x,fx,'o');
hold on;
plot(xx,ffx,'r');
hold on;
fx=1./(1+25*(xx.^2));
plot(xx,fx,'k');
hold on;
difference=ffx-fx;
plot(xx,difference,'g');
title('三次样条插值----三转角法');
legend('节点','插值结果','原函数','误差');
grid on;