mathanalysis303.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))+h*diag(ones(N-1,1),1)+h*diag(ones(N-1,1),-1);
H=H(2:N,2:N);
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(1)=D(1)-h*ddfx(1)/6;
D(N-1)=D(N-1)-h*ddfx(N-1)/6;
D=6*D;
M=mathanalysis202(H,D);
M=[ddfx(1)
    M
    ddfx(N-1)];
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;