上交Newton拉格朗日插值MATLAB大作业编程.docx
- 文档编号:26946209
- 上传时间:2023-06-24
- 格式:DOCX
- 页数:13
- 大小:61.09KB
上交Newton拉格朗日插值MATLAB大作业编程.docx
《上交Newton拉格朗日插值MATLAB大作业编程.docx》由会员分享,可在线阅读,更多相关《上交Newton拉格朗日插值MATLAB大作业编程.docx(13页珍藏版)》请在冰豆网上搜索。
上交Newton拉格朗日插值MATLAB大作业编程
作业报告
使用MATLAB编程
2.1
Newton插值:
clear
clc
x=[0.20.40.60.81];
n=length(x);
y=zeros(n);
y(:
1)=[0.980.920.810.640.38];
a=x
(1);
b=x(end);
forj=1:
n-1
fork=j:
n-1
temp=y(k+1,j)-y(k,j);
y(k+1,j+1)=temp/(x(k+1)-x(k+1-j));
end
c(j)=y(j,j);
c(j+1)=y(j+1,j+1);
end
p=c
(1);
q=1;
symsX
fori=2:
n
q=q*(X-x(i-1));
p=p+c(i)*q;
end
p=simple(p)
fori=1:
301
t(i)=a+(b-a)/300*(i-1);
Nn(i)=eval(subs(p,'X','t(i)'));
end
plot(t,Nn,'b')
holdon
gridon
legend('Newton')
title('Newton')
xlabel('x')
ylabel('f(x)')
图像:
三次样条插值:
M文件spline33.m
function[yybcd]=spline3(x,y,xx,flag,vl,vr)
iflength(x)~=length(y)
error('输入数据应成对');
end
n=length(x);
a=zeros(n-1,1);
b=a;d=a;dx=a;dy=a;
A=zeros(n);B=zeros(n,1);
fori=1:
n-1
a(i)=y(i);
dx(i)=x(i+1)-x(i);
dy(i)=y(i+1)-y(i);
end
fori=2:
n-1
A(i,i-1)=dx(i-1);
A(i,i)=2*(dx(i-1)+dx(i));
A(i,i+1)=dx(i);
B(i,1)=3*(dy(i)/dx(i)-dy(i-1)/dx(i-1));
end
ifflag==0
A(1,1)=1;
A(n,n)=1;
end
%---------------------------------%
%端点一阶导数条件%
ifflag==1
A(1,1)=2*dx
(1);A(1,2)=dx
(1);
A(n,n-1)=dx(n-1);A(n,n)=2*dx(n-1);
B(1,1)=3*(dy
(1)/dx
(1)-vl);
B(n,1)=3*(vr-dy(n-1)/dx(n-1));
end
%---------------%
%端点二阶导数条件%
ifflag==2
A(1,1)=2;A(n,n)=2;
B(1,1)=vl;B(n,1)=vr;
end
%---------------%
c=A\B;
fori=1:
n-1
d(i)=(c(i+1)-c(i))/(3*dx(i));
b(i)=dy(i)/dx(i)-dx(i)*(2*c(i)+c(i+1))/3;
end
[mmnn]=size(xx);
yy=zeros(mm,nn);
fori=1:
mm*nn
forii=1:
n-1
ifxx(i)>=x(ii)&xx(i) j=ii; break; elseifxx(i)==x(n) j=n-1; end end yy(i)=a(j)+b(j)*(xx(i)-x(j))+c(j)*(xx(i)-x(j))^2+d(j)*(xx(i)-x(j))^3; end end 主程序: clear; clc; formatshortg; x1=[0.20.40.60.81.0]'; y1=[0.980.920.810.640.38]'; u1=0;un=0; xx1=[x1 (1): 0.001: x1(end)]'; [yy1b1c1d1]=spline33(x1,y1,xx1,0,u1,un); fprintf('\t\tb1\t\tc1\t\td1\n'); disp([b1c1(1: end-1,1)d1]); plot(x1,y1,'bo',xx1,yy1,'--k'); legend('spline') title('spline') gridon; 2.3 (1)8次拉格朗日插值 M文件language1.m functionf=language(x,y,x0) symstl; if(length(x)==length(y)) n=length(x); else disp('x和y的维数不相等! ') return; end h=sym(0); fori=1: n l=sym(y(i)); forj=1: i-1 l=l*(t-x(j))/(x(i)-x(j)); end; forj=i+1: n l=l*(t-x(j))/(x(i)-x(j)); end; h=h+l; end simplify(h); if(nargin==3) f=subs(h,t,x0); else f=collect(h); f=vpa(f,6); end 主程序: x=[01491625364964]; y=[012345678]; f=language1(x,y) plot(x,y) legend('language') title('language') gridon; (2)三次样条差值: M文件spline33.m function[yybcd]=spline3(x,y,xx,flag,vl,vr) iflength(x)~=length(y) error('输入数据应成对'); end n=length(x); a=zeros(n-1,1); b=a;d=a;dx=a;dy=a; A=zeros(n);B=zeros(n,1); fori=1: n-1 a(i)=y(i); dx(i)=x(i+1)-x(i); dy(i)=y(i+1)-y(i); end fori=2: n-1 A(i,i-1)=dx(i-1); A(i,i)=2*(dx(i-1)+dx(i)); A(i,i+1)=dx(i); B(i,1)=3*(dy(i)/dx(i)-dy(i-1)/dx(i-1)); end ifflag==0 A(1,1)=1; A(n,n)=1; end %---------------------------------% %端点一阶导数条件% ifflag==1 A(1,1)=2*dx (1);A(1,2)=dx (1); A(n,n-1)=dx(n-1);A(n,n)=2*dx(n-1); B(1,1)=3*(dy (1)/dx (1)-vl); B(n,1)=3*(vr-dy(n-1)/dx(n-1)); end %---------------% %端点二阶导数条件% ifflag==2 A(1,1)=2;A(n,n)=2; B(1,1)=vl;B(n,1)=vr; end %---------------% c=A\B; fori=1: n-1 d(i)=(c(i+1)-c(i))/(3*dx(i)); b(i)=dy(i)/dx(i)-dx(i)*(2*c(i)+c(i+1))/3; end [mmnn]=size(xx); yy=zeros(mm,nn); fori=1: mm*nn forii=1: n-1 ifxx(i)>=x(ii)&xx(i) j=ii; break; elseifxx(i)==x(n) j=n-1; end end yy(i)=a(j)+b(j)*(xx(i)-x(j))+c(j)*(xx(i)-x(j))^2+d(j)*(xx(i)-x(j))^3; end end 主程序: clear; clc; formatshortg; x1=[01491625364964]'; y1=[012345678]'; u1=0;un=0; xx1=[x1 (1): 0.001: x1(end)]'; [yy1b1c1d1]=spline33(x1,y1,xx1,1,u1,un); fprintf('\t\tb1\t\tc1\t\td1\n'); disp([b1c1(1: end-1,1)d1]); plot(x1,y1,'bo',xx1,yy1,'--k'); legend('spline') title('spline') gridon; 4-1 (1)M文件function.m: Functiony=fun(x); y=sqrt(x)*log(x); 主程序 Clear;clc; h=0.001; n=1/h;t=0;s1=0;s2=0; fori=1: n-1 t=t+function(i*h); end T=h/2*(0+2*t+f (1)); T=vpa(T,7) Fori=0: n-1 s1=s1+function(h/2+i*h); end fori=1: n-1 s2=s2+function(i*h); end S=h/6*(0+4*s1+2*s2+f (1));S=vpa(S,7) 综上,得解。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 上交 Newton 拉格朗日插值 MATLAB 作业 编程