else
disp('y和y的导数的维数不相等！'); return; end else
disp('x和y的维数不相等！'); return; end
for i=1:n h = 1.0; a = 0.0; for j=1:n if( j ~= i)
h = h*(t-x(j))^2/((x(i)-x(j))^2); a = a + 1/(x(i)-x(j)); end end
f = f + h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i));
if(i==n)
if(nargin == 4)
f = subs(f,'t',x0); else
f = vpa(f,6); end end End
10.Language
function f = Language(x,y,x0) syms t;
if(length(x) == length(y)) n = length(x); else
disp('x和y的维数不相等！'); return;
end
f = 0.0; for(i = 1:n) l = y(i);
%检错 for(j = 1:i-1)
l = l*(t-x(j))/(x(i)-x(j)); end;
for(j = i+1:n)
l = l*(t-x(j))/(x(i)-x(j)); %计算拉格朗日基函数 end;
f = f + l; %计算拉格朗日插值函数 simplify(f); %化简
if(i==n)
if(nargin == 3)
f = subs(f,'t',x0); %计算插值点的函数值 else
f = collect(f); %将插值多项式展开
f = vpa(f,6); %将插值多项式的系数化成6位精度的小数 end end end
11.Neville
function f = Neville(x,y,x0) syms t;
if(length(x) == length(y)) n = length(x); else
disp('x和y的维数不相等！'); return; end
y1(1:n) = t; for(i=1:n-1) for(j=i+1:n) if(j==2)
y1(j) = y(j)+(y(j)-y(j-1))/((t-x(j-i))/(t-x(j)))*(1-(y(j)-y(j-1))/y(j)); else
y1(j) = y(j)+(y(j)-y(j-1))/((t-x(j-i))/(t-x(j)))*(1-(y(j)-y(j-1))/(y(j)-y(j-2))); end
end y = y1; if(i==n-1)
if(nargin == 3)
f = subs(y(n-1),'t',x0);
else
f = vpa(y(n-1),6); end end end
12.Newton
function f = Newton(x,y,x0) syms t;
if(length(x) == length(y)) n = length(x); c(1:n) = 0.0; else
disp('x和y的维数不相等！'); return; end
f = y(1); y1 = 0; l = 1;
for(i=1:n-1) for(j=i+1:n)
y1(j) = (y(j)-y(i))/(x(j)-x(i)); end
c(i) = y1(i+1); l = l*(t-x(i)); f = f + c(i)*l; simplify(f); y = y1;
if(i==n-1)
if(nargin == 3)
f = subs(f,'t',x0); else
f = collect(f); f = vpa(f, 6); end end end
%将插值多项式展开
13.Newtonback
function f = Newtonback(x,y,x0) syms t;
if(length(x) == length(y)) n = length(x); c(1:n) = 0.0; else
disp('x和y的维数不相等！'); return; end
f = y(n); y1 = 0;
xx =linspace(x(1),x(n),(x(2)-x(1))); if(xx ~= x)
disp('节点之间不是等距的！'); return; end
for(i=1:n-1) for(j=i+1:n)
y1(j) = y(j)-y(j-1); end
c(i) = y1(n); l = t;
for(k=1:i-1) l = l*(t+k); end;
f = f + c(i)*l/factorial(i); simplify(f); y = y1;
if(i==n-1)
if(nargin == 3)
f = subs(f,'t',(x0-x(n))/(x(2)-x(1))); else
f = collect(f); f = vpa(f, 6); end end
end
14.Newtonforward
function f = Newtonforward(x,y,x0) syms t;
if(length(x) == length(y)) n = length(x); c(1:n) = 0.0; else
disp('x和y的维数不相等！'); return; end
f = y(1); y1 = 0;
xx =linspace(x(1),x(n),(x(2)-x(1))); if(xx ~= x)
disp('节点之间不是等距的！'); return; end
for(i=1:n-1) for(j=1:n-i)
y1(j) = y(j+1)-y(j); end
c(i) = y1(1); l = t;
for(k=1:i-1) l = l*(t-k); end;
f = f + c(i)*l/factorial(i); simplify(f); y = y1;
if(i==n-1)
if(nargin == 3)
f = subs(f,'t',(x0-x(1))/(x(2)-x(1))); else
f = collect(f); f = vpa(f, 6); end end