数值分析的lagrange插值为例讲讲如何在北太天元软件上提速
%计算x0处的插值函数L(x)的值
% @param [in] x: 插值点的x坐标
% @param [in] y: 插值点的y坐标
% @param [in] x0: 要计算x0处的函数值L(x0), 现在的x0是1x1 double, 后续
% 希望能改成x0 可以是 1xn 的矩阵
% @param [out] z: 返回值z=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 是一个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 之间的分点，均匀剖分，分成100等分（产生101个点)
% f 是在这些剖分点上的值 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 表示矩阵的逐个分量做的平方
%% 计算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份取得的插值点
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('N=%d\n',N);
fprintf('Max Error of grid (1):%f\n',d1);
fprintf('Max Error of grid (2):%f\n',d2);
end
%%再次计算N=10时的插值函数，并且画图
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');
