|
标签（空格分隔）： matlab 插值 算法 数值 分析 lagrange Newton Hermite spline
matlab lagrange插值，Newton插值，Hermite插值，spline插值
1 lagrange插值
1.1 测试lagrange插值函数lagrangeInsertValue()与lagrangeInsertValueWithError()
clc,clear,clf
[email protected](x)(sqrt(x));
% x=[4 9 6.25 40]
x=[100 121 81 144]
y=f(x)
% xs=[7 8 100]
xs=[105]
% x=[-1 0 .5 1 ]
% y=[3 -.5 0 1]
% syms xs
%绘制插值函数图
hold on
plot(x,y,'o')
min_=min(x);max_=max(x);
xs_pic=min_:(max_-min_)/50:max_;
ys_pic=lagrangeInsertValue(x,y,xs_pic);
plot(xs_pic,ys_pic,'-.')
ys_pic=f(xs_pic);
plot(xs_pic,ys_pic)
xs_pic=[];ys_pic=[];
hold off
%%%%%%%%%%%%计算插值、误差
% ys=lagrangeInsertValue(x,y,xs)
%增加的用于求误差的点
xadd=prod(xs)^(1/length(xs));xadd=xadd+.1*max(xs);
yadd=f(xadd);
[ys,rs]=lagrangeInsertValueWithError(x,y,xs,xadd,yadd)
xadd,yadd
rs_true=f(xs)-ys
得到函数曲线与插值曲线如下：
1.2 lagrangeInsertValue()
lagrange插值函数，不包含误差计算功能
function ys=lagrangeInsertValue(x,y,xs)
%【程序特点】：考虑将内存和运算量同时达到极小化
%ys为插值xs的函数值
%rs为对应的插值误差
%时间：2013/12/08 姓名：邓能财
%%% 标记为显示行（可去）%#符号
len=length(x);
vlen=length(xs);
% %计算xi-xj形成的矩阵
% dx=zeros(len);
% for i=1:len
% for j=i+1:len
% dx(i,j)=x(i)-x(j);
% end
% end
% dx=dx-dx'; %xj-xi=-(xi-xj)
% for i=1:len
% dx(i,i)=1;
% end
%dx
%%%%%%%%%%%%计算w(xi)
wxi=ones(1,len);
for i=1:len
for j=1:len
if j~=i
wxi(i)=wxi(i)*(x(i)-x(j));
end
end
end
wxi
%%%%%%%%%%%%计算函数值
%用【卷计算的方式节约到vlen个double内存】
cmMultidx=zeros(1,vlen);
for j=1:vlen
cmMultidx(j)=prod(xs(j)-x);
end
cmMultidx%#
ys=zeros(1,vlen);
sum_l=ys;%#
for i=1:len
l=cmMultidx./((xs-x(i))*wxi(i))
sum_l=sum_l+l;%#
ys=ys+l.*y(i);
end
sum_l%#
end
1.3 lagrangeInsertValueWithError()
lagrange插值函数，包含误差计算功能
function [ys,rs]=lagrangeInsertValueWithError(x,y,xs,xadd,yadd)
%lagrange插值，包含【事后误差估计】版
%时间：2013/12/08 姓名：邓能财
ys=lagrangeInsertValue(x,y,xs);
rs=ys-lagrangeInsertValue([x(2:end),xadd],[y(2:end),yadd],xs);
end
2 Newton插值
2.1 测试Newton插值函数newtonInsertValue()
clc,clear,clf
% x=[4 9 6.25]
% y=[2 3 2.5]
% xs=[7 8 8.5]
[email protected](x)(x.^5+x.^3+1);
x=[-1 -.8 0 .5 1]
y=f(x)
%绘制插值函数图
hold on
plot(x,y,'o')
min_=min(x);max_=max(x);
xs_pic=min_:(max_-min_)/50:max_;
ys_pic=newtonInsertValue(x,y,xs_pic);
plot(xs_pic,ys_pic,'-.')
ys_pic=f(xs_pic);
plot(xs_pic,ys_pic)
xs_pic=[];ys_pic=[];
hold off
xs=[-.4 .7 4]
[ys,rs]=newtonInsertValue(x,y,xs)
xs
得到函数曲线与插值曲线如下：
2.2 newtonInsertValue()
function [ys,rs]=newtonInsertValue(x,y,xs)
%【程序特点】：考虑将内存和运算量同时达到极小化
%ys为插值xs的函数值
%rs为对应的插值误差
%时间：2013/12/08 姓名：邓能财
%%% 标记为显示行（可去）%#符号
len=length(x);
vlen=length(xs);
%%%%%%%%计算各阶差商，在计算差商中只记录：
%1 插值多项式系数的差商df1
%2 用于误差计算的差商dfn
%3 计算过程中间的i阶差商dfi
df1=zeros(1,len); dfn=df1; dfi=df1;
dfi=y;
dfi_show=dfi(1:len)%#
df1(1)=dfi(1)
dfn(1)=dfi(len)
for i=1:len-1 %差商阶数
for j=1:len-i %列标
dfi(j)=(dfi(j+1)-dfi(j))/(x(j+i)-x(j));
end
dfi_show=dfi(1:len-i)%#
df1(i+1)=dfi(1)
dfn(i+1)=dfi(len-i)
end
dfi=[]; %释放内存
%%%%%%%%%计算插值多项式的值
%用【卷计算的方式节约到vlen个double内存】
%设置卷变量：wxi=(xs-x1)*...*(xs-xi)
wxi=ones(1,vlen);
ys=df1(1)*wxi;
for i=2:len
wxi=wxi.*(xs-x(i-1));
ys=ys+df1(i)*wxi;
end
df1=[]; %释放内存
%%%%%%%%%计算误差项
%计算包含xs的差商df_nx
dfn_x=ys; %df_nx
for i=1:len
dfn_x=(dfn_x-dfn(i))./(xs-x(len+1-i));
end
if length(dfn_x)<10, dfn_x,end%#
wxi=wxi.*(xs-x(len));
rs=dfn_x.*wxi;
end
3 Hermite插值
3.1 测试Hermite插值函数hermiteInsertValue()
clc,clear,clf
%%%%%%%%e1
% [email protected](x)(sqrt(x));
% [email protected](x)(.5*x.^(-.5));
% x=[4 9 6.25 40]
% xs=[7 8 100]
%%%%%%%%%e2
[email protected](x)(sqrt(x).*sin(x));
[email protected](x)(.5*x.^(-.5).*sin(x)+sqrt(x).*cos(x));
x=[4 9 16.5]
y=f(x)
dy=df(x)
xs=[7 8 20]
%绘制插值函数图
hold on
plot(x,y,'o')
min_=min(x);max_=max(x);
h_pic=(max_-min_)/50;
xs_pic=min_-1:h_pic:max_+1;
ys_pic=hermiteInsertValue(x,y,dy,xs_pic);
plot(xs_pic,ys_pic,'-.')
ys_pic=f(xs_pic);
plot(xs_pic,ys_pic)
xs_pic=[];ys_pic=[];
hold off
%%%%%%
disp('%%%%%%%%%%%%计算插值、误差')
% [ys,rs]=
ys=hermiteInsertValue(x,y,dy,xs)
得到函数曲线与插值曲线如下：
3.2 hermiteInsertValue()
function ys=hermiteInsertValue(x,y,dy,xs)
%【程序特点】：考虑将内存和运算量同时达到极小化
%ys为插值xs的函数值
%rs为对应的插值误差
%时间：2013/12/08 姓名：邓能财
%%% 标记为显示行（可去）%#符号
max_show_numb=10;
len=length(x);
vlen=length(xs);
%%%%%%%%%%%%计算w(xi)、sum_dx(xi)
wxi=ones(1,len);
sum_dxi=zeros(1,len);
for i=1:len
for j=1:len
if j~=i
wxi(i)=wxi(i)*(x(i)-x(j));
sum_dxi(i)=sum_dxi(i)+1/(x(i)-x(j));
end
end
end
if len<max_show_numb, wxi,sum_dxi, end%#
%%%%%%%%%%%%计算函数值
%用【卷计算的方式节约到vlen个double内存】
cmMultidx=zeros(1,vlen);
for j=1:vlen
cmMultidx(j)=prod(xs(j)-x);
end
if vlen<max_show_numb, cmMultidx, end%#
ys=zeros(1,vlen);
sum_l=ys;%#
for i=1:len
l=cmMultidx./((xs-x(i))*wxi(i));
sum_l=sum_l+l;%#
gx_f=(1-2*(xs-x(i))*sum_dxi(i)).*l.*l;
hx_f_=(xs-x(i)).*l.*l;
ys=ys+gx_f.*y(i)+hx_f_.*dy(i);
end
if vlen<max_show_numb, sum_l, end%#
end
4 spline插值
4.1 测试spline插值函数splineInterpolation()
clc,clear;
clf
% [email protected](x)(x.^5+x.^3+1);
% [email protected](x)(sqrt(x).*sin(x));
[email protected](x)(log(x).*sin(x));
% x=[-1 -.8 0 .5 1]
% x=[.25 .30 .39 .45 .53]
% y=[.50 .54 .62 .67 .73]
x=[1 3 5.2 7.7 9]
y=f(x)
% % 绘制插值函数图
hold on
plot(x,y,'o')
min_=min(x);max_=max(x);
h_pic=(max_-min_)/50; d_pic=6;
xs_pic=min_-d_pic*h_pic:h_pic:max_+d_pic*h_pic;
ys_pic=splineInterpolation(x,y,xs_pic);
plot(xs_pic,ys_pic,'-.')
ys_pic=f(xs_pic);
plot(xs_pic,ys_pic)
xs_pic=[];ys_pic=[];
hold off
%%%%%%
xs=[-.4 .7 4]
m0=1,mn=.68;
ys=splineInterpolation(x,y,xs,m0,mn)
xs
得到函数曲线与插值曲线如下：
4.2 splineInterpolation()
function [ys,rs]=splineInterpolation(x,y,xs,mstart,mend)
%【程序特点】：考虑将内存和运算量同时达到极小化
%ys为插值xs的函数值
%rs为对应的插值误差
%时间：2013/12/08 姓名：邓能财
%%% 标记为显示行（可去）%#符号
%%%%%%%%%%%%默认端点条件为：两端导数为0
if nargin<4, mstart=0;mend=0; end
len=length(x);
vlen=length(xs);
%%%%%%%%%%%%将x从小到大排序，y也随之排序
[x,index]=sort(x); y=y(index); index=[];
%%%%%%%%%%%%计算一阶导数m
%计算线性方程的系数
lamda=zeros(1,len-2); mu=lamda; c=lamda;
i=2;
df_memory=(y(i)-y(i-1))/(x(i)-x(i-1)); df=df_memory;
hi_memory=x(i)-x(i-1); hi=hi_memory;
for im=1:len-2
i=im+1;
hi_memory=hi; hi=x(i+1)-x(i);
df_memory=df; df=(y(i+1)-y(i))/hi;
lamda(im)=hi/(hi_memory+hi); mu(im)=1-lamda(im);
c(im)=3*(lamda(im)*df_memory+mu(im)*df);
end
lamda,mu,c%#
%求解线性方程，计算m
c(1)=c(1)-lamda(1)*mstart; c(end)=c(end)-mu(end)*mend;
m=zeros(1,len); m(1)=mstart; m(end)=mend;
m(2:end-1)=bandMatrixSolve(2*ones(1,len-2),mu(1:end-1),lamda(2:end),c);
m%#
%%%%%%%%%%%%计算函数值
%用【卷计算的方式节约到vlen个double内存】
%x已排序
%将xs排序,再查找xs(i)所在的区间[xi-1,xi]
[xs,index]=sort(xs);
ys=zeros(1,vlen);
j=1;
for i=1:vlen
%从当前位置（非起始）开始查找xs(i)所在的区间
while j+1<len && xs(i)>x(j+1), j=j+1; end
%计算ys
h=x(j+1)-x(j); h3_=1/(h*h*h);
dx1=xs(i)-x(j); dx2=dx1*dx1;
dx1_=xs(i)-x(j+1); dx2_=dx1_*dx1_;
ys(index(i))=h3_*h*(dx2_*y(j)+dx2*y(j+1)...
+dx1*dx2_*m(j)+dx1_*dx2*m(j+1))...
+h3_*(2*dx1*dx2_*y(j)...
-2*dx1_*dx2*y(j+1));
end
rs=0;
end
4.3 bandMatrixSolve()
function x=bandMatrixSolve(a,b,c,e)
%【程序特点】：考虑将内存和运算量同时达到极小化
%时间：2013/12/08 姓名：邓能财
%带状矩阵方程组求解
%矩阵A由中a(dim维)、b、c，构成
%等式右边的向量为：e
%算法：追赶法
dim=length(a);
assert( dim==length(b)+1 &&dim==length(c)+1 &&dim==length(e),...
['Argument input error: ',...
'Matrix dimensions must agree.'])
%%%%%%%%%%%%%%%%%%%LU分解
alfa=zeros(1,dim); beta=zeros(1,dim-1);
%初始值
alfa(1)=a(1);
beta(1)=b(1)/alfa(1);
for i=2:dim-1
alfa(i)=a(i)-c(i-1)*beta(i-1);
beta(i)=b(i)/alfa(i);
end
i=i+1;
alfa(i)=a(i)-c(i-1)*beta(i-1);
a=[]; b=[];%释放内存
%%%%%%%%%%%%%%%%%%%解Ly=e
y=zeros(1,dim);
%初始值
y(1)=e(1)/alfa(1);
for i=2:dim
y(i)=(e(i)-c(i-1)*y(i-1))/alfa(i);
end
alfa=[]; c=[];%释放内存
%%%%%%%%%%%%%%%%%%%解Ux=y
x=zeros(1,dim);
%初始值
x(dim)=y(dim);
for i=dim-1:-1:1
x(i)=y(i)-beta(i)*x(i+1);
end
beta=[]; y=[];%释放内存
end
end
|
客服电话
APP下载
官方微信
请发表评论