数值分析拟合MATLAB代码
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第二章
第一题
MATLAB代码：
%用spline作图
xi=[0.20.40.60.81.0];
yi=[0.980.920.810.640.38];
x=-1:0.01:2;
y1=Newton3(xi,yi,x);%源代码见m文件y2=spline(xi,yi,x);
plot(xi,yi,'o',x,y1,'r-',x,y2,'k')
%用csape作图
xi=[0.20.40.60.81.0];
yi=[0.980.920.810.640.38];
x=-1:0.01:2;
y1=Newton3(xi,yi,x);%源代码见m文件pp=csape(xi,yi,'variational');
y2=ppval(pp,x);
plot(xi,yi,'o',x,y1,'r-',x,y2,'k')
%Newton3
function f=Newton3(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
第二题
MATLAB代码：
y1=zeros(1,11);
x1=-1:0.2:1;
for i=1:11,
p=1./(1+25.*x1(i).^2);
y1(i)=p;
end
x=-1:0.01:1;
ym=language(x1,y1,x);
yn=spline(x1,y1,x);
y2=zeros(1,21);
x2=-1:0.1:1;
for i=1:21,
p=1./(1+25.*x2(i).^2);
y2(i)=p;
end
x=-1:0.01:1;
yi=language(x2,y2,x);
yj=spline(x2,y2,x);
figure(1);
plot(x1,y1,'o',x,ym,'r-',x,yn,'k-'); figure(2);
plot(x2,y2,'o',x,yi,'r-',x,yj,'k-')
%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);
end
end
end
第三题
MATLAB代码：
用spline进行插值
x1=[01491625364964]; y1=[012345678];
x=0:0.1:64;
ym=language(x1,y1,x);
yn=spline(x1,y1,x);
figure(1);
plot(x1,y1,'o',x,ym,'k-',x,yn,'r-');
用scape用第一边界条件进行插值
x1=[01491625364964];
y1=[012345678];
x=0:0.1:64;
ym=language(x1,y1,x);%原函数见上面pp=csape(x1,y1,'complete',[0.2,-1]);
yn=ppval(pp,x);
figure(1);
plot(x1,y1,'o',x,ym,'k-',x,yn,'r-');
结论：
[064]三次样条插值精确
[01]多项式插值精确
第三章
第一题
MATLAB代码：
y=zeros(1,11);
x=-1:0.2:1;
for i=1:11,
p=1./(1+25.*x(i).^2);
y(i)=p;
end
A=polyfit(x,y,10);
y1=poly2str(A,'x');
p1=-1:0.02:1;
u1=polyval(A,p1);
plot(p1,u1,x,y,'o')