%計(jì)算x0處的插值函數(shù)L(x)的值

% @param [in] x: 插值點(diǎn)的x坐標(biāo)

% @param [in] y: 插值點(diǎn)的y坐標(biāo)

% @param [in] x0: 要計(jì)算x0處的函數(shù)值L(x0), 現(xiàn)在的x0是1x1 double, 后續(xù)

% 希望能改成x0 可以是 1xn 的矩陣

% @param [out] z: 返回值z(mì)=L(x0)

function z=lagrint(x,y,x0)

n=length(x);

l=ones(1,n);

omega_x = x0 - x;

for i=1:n

/*

for j=1:n

if j~=i

l(i)=l(i)*(x0-x(j))/(x(i)-x(j));

end

end

*/

x_ij = x(i)- x;

if i == 1

l(i) = prod(omega_x(i+1:n)) / prod(x_ij(i+1:n));

elseif i == n

l(i) = prod(omega_x(1:i-1)) / prod(x_ij(1:i-1));

else

l(i) = prod(omega_x(1:i-1)) / prod(x_ij(1:i-1));

l(i) = l(i)*prod(omega_x(i+1:n)) / prod(x_ij(i+1:n));

end

end

z=y*l';

end

% n1 是一個(gè)1x4的矩陣，取值分別是 5,10,20,40;

%n1=ones(1,4);

%n1(1)=5;

%n1(2)=10;

%n1(3)=20;

%n1(4)=40;

%在北太天元上更加高效的寫法

n1 = [5,10,20,40];

%

% y = -5,5 之間的分點(diǎn)，均勻剖分，分成100等分（產(chǎn)生101個(gè)點(diǎn))

% f 是在這些剖分點(diǎn)上的值 f(t) = 1/(1+t^2) t= y(1), ..., y(101)

/*

y=ones(1,101);

f=ones(1,101);

for i=1:101

???y(i)=-5+(i-1)/10;

???f(i)=1/(1+y(i)^2);

end

*/

y = linspace(-5,5,101);

f = 1 / ( 1 +?real(y.^2) ); %?.^2 表示矩陣的逐個(gè)分量做的平方

%% 計(jì)算N=5,10,20,40 的插值誤差

p1=ones(1,101);

p2=ones(1,101);

for i=1:4

???N=n1(i);

???x1=ones(1,N+1);

???x2=ones(1,N+1);

???f1=ones(1,N+1);

???f2=ones(1,N+1);

???for j=1:(N+1)

???????x1(j)=5-10*(j-1)/N; %[-5,5]平均分成N份取得的插值點(diǎn)

???????x2(j)=-5*cos(pi*(2*j-1)/(2*N+2)); %切比雪夫插值點(diǎn)

???????f1(j)=1/(1+x1(j)^2); %均分插值點(diǎn)上的函數(shù)值

???????f2(j)=1/(1+x2(j)^2); %切比雪夫插值點(diǎn)上的函數(shù)值

???end

???d1=0;

???d2=0;

???for j=1:101

???????p1(j)=lagrint(x1,f1,y(j));

???????p2(j)=lagrint(x2,f2,y(j));

???????if d1<abs(f(j)-p1(j))

???????????d1=abs(f(j)-p1(j));

???????end

???????if d2<abs(f(j)-p2(j))

???????????d2=abs(f(j)-p2(j));

???????end

???end

???fprintf('N=%d\n',N);

???fprintf('Max Error of grid (1):%f\n',d1);

???fprintf('Max Error of grid (2):%f\n',d2);

end

%%再次計(jì)算N=10時(shí)的插值函數(shù)，并且畫圖

N=10;

x1=ones(1,N+1);

x2=ones(1,N+1);

f1=ones(1,N+1);

f2=ones(1,N+1);

for j=1:(N+1)

???x1(j)=5-10*(j-1)/N;

???x2(j)=-5*cos(pi*(2*j-1)/(2*N+2));

???f1(j)=1/(1+x1(j)^2);

???f2(j)=1/(1+x2(j)^2);

end

d1=0;

d2=0;

for j=1:101

???p1(j)=lagrint(x1,f1,y(j));

???p2(j)=lagrint(x2,f2,y(j));

???if d1<abs(f(j)-p1(j))

???????d1=abs(f(j)-p1(j));

???end

???if d2<abs(f(j)-p2(j))

???????d2=abs(f(j)-p2(j));

???end

end

fprintf('The graph of N=10(r:origin,b:grid 1,g:grid 2):\n');

plot(y,f,'r',y,p1,'b',y,p2,'g');