数值分析作业.docx
- 文档编号:16726020
- 上传时间:2023-07-16
- 格式:DOCX
- 页数:32
- 大小:373.78KB
数值分析作业.docx
《数值分析作业.docx》由会员分享,可在线阅读,更多相关《数值分析作业.docx(32页珍藏版)》请在冰点文库上搜索。
数值分析作业
数值分析实验作业
指导老师:
院系:
姓名:
学号:
第二章
用matlab编写的程序如下:
functiont_charpt2
%数值实验二:
含“实验2.1:
多项式插值的震荡现象”和“实验2.2:
样条插值的收敛”
%输入:
实验选择,函数式选择,插值结点数
%输出:
拟合函数及原函数的图形
result=inputdlg({'请选择实验,若选2.1,请输入1,否则输入2:
'},'charpt2',1,{'1'});
Nb=str2num(char(result));
if(Nb~=1)&(Nb~=2)errordlg('实验选择错误!
');return;end
promps={'请选择实验函数,若选f(x),请输入f,若选h(x),请输入h,若选g(x),请输入g:
'};
titles='charpt2';
result=inputdlg(promps,'charpt2',1,{'f'});
Nb_f=char(result);
if(Nb_f~='f'&Nb_f~='h'&Nb_f~='g')errordlg('实验函数选择错误!
');return;end
result=inputdlg({'请输入插值结点数N:
'},'charpt2',1,{'10'});
Nd=str2num(char(result));
if(Nd<1)errordlg('结点输入错误!
');return;end
switchNb_f
case'f'
f=inline('1./(1+25*x.^2)');a=-1;b=1;
case'h'
f=inline('x./(1+x.^4)');a=-5;b=5;
case'g'
f=inline('atan(x)');a=-5;b=5;
end
if(Nb==1)
x0=linspace(a,b,Nd+1);y0=feval(f,x0);
x=a:
0.1:
b;y=Lagrange(x0,y0,x);
fplot(f,[ab],'co');
holdon;
plot(x,y,'b--');
xlabel('x');ylabel('y=f(x)oandy=Ln(x)--');
elseif(Nb==2)
x0=linspace(a,b,Nd+1);y0=feval(f,x0);
x=a:
0.1:
b;
cs=spline(x0,y0);y=ppval(cs,x);
plot(x0,y0,'o');holdon;plot(x,y,'k-');
xlabel('x');ylabel('y=f(x)oandy=Spline(x)-');
end
%
functiony=Lagrange(x0,y0,x)
n=length(x0);m=length(x);
fori=1:
m;
z=x(i);
s=0.0;
fork=1:
n
p=1.0;
forj=1:
n
if(j~=k)
p=p*(z-x0(j))/(x0(k)-x0(j));
end
end
s=s+p*y0(k);
end
y(i)=s;
end
若选择实验2.1,实验函数为f,插值结点数为6,则结果为:
若选择实验2.1,实验函数为f,插值结点数为20,则结果为:
若选择实验2.1,实验函数为f,插值结点数为27,则结果为:
选择其它的函数重复上述的实验,在这里我选择的是h函数,具体结果如下:
若选择实验2.1,实验函数为h,插值结点数为6,则结果为:
若选择实验2.1,实验函数为h,插值结点数为20,则结果为:
若选择实验2.1,实验函数为h,插值结点数为27,则结果为:
实验2.2
①输入2,实验函数为f,插值节点数为20,则结果为:
②输入2,实验函数为h,插值节点数为20,则结果为:
分析:
1.对于n次Lagrange插值多项式Pn(x)近似逼近函数f(x)时,Pn(x)的次数并非越高,越逼近目标函数f(x)。
在实验中选取一个较小数6,可以看到Pn(x)图像与f(x)吻合情况不是很好,当选取n=20时,Pn(x)与f(x)在区间[-1,1]上的吻合情况较好;但是当选取较大数27,Pn(x)开始发生振荡,吻合情况较差,也即“龙格”现象。
在选取的函数h(x)不会发生类似的现象,精度会随着次数的增加而增加,不会出现所谓的“龙格”现象。
2.从样条插值的实验结果可以看到:
样条插值的逼近情况比Lagrange插值要好,并且Lagrange插值在区间[-5,5]上端点处会有一定的波动,但是样条插值结果与精确值吻合较好。
第三章实验3
用matlab编写的程序如下:
functioncharpt3
%数值实验三:
含“实验3.1”和“实验3.2”
%子函数调用:
dlsa
%输入:
实验选择
%输出:
原函数及求得的相应插值多项式的函数的图像以及参数alph和误差r
result=inputdlg({'请选择实验,若选3.1,请输入1,否则输入2:
'},'charpt3',1,{'1'});
Nb=str2num(char(result));
if(Nb~=1)&(Nb~=2)errordlg('实验选择错误!
');return;end
x0=-1:
0.5:
2;
y0=[-4.447-0.4520.5510.048-0.4470.5494.552];
if(Nb==1)
n=3;%n为拟合阶次
alph=polyfit(x0,y0,n);
y=polyval(alph,x0);
r=(y0-y)*(y0-y)';
x=-1:
0.01:
2;
y=polyval(alph,x);
plot(x,y,'k--');
xlabel('x');ylabel('y0*andployfit.y--');
holdon
plot(x0,y0,'*')
gridon;
else
result=inputdlg({'请输入权向量w:
'},'charpt3',1,{'[1111111]'});
w=str2num(char(result));
n=3;
[a,b,c,alph,r]=dlsa(x0,y0,w,n);
end
disp(['平方误差:
',num2str(r)])
disp(['参数alph:
',num2str(alph)])
function[a,b,c,alph,r]=dlsa(x,y,w,n)
%功能:
用正交化方法对离散数据作多项式最小二乘拟合。
%输入:
m+1个离散点(x,y,w),x,y,w分别用行向量给出。
%拟合多项式的次数n,0 %平方误差r=(y-s,y-s),并作离散点列和拟合曲线的图形 m=length(x)-1; if(n<1|n>=m)errordlg('错误: n<1或者n>=m! ');return;end %求三项递推公式的参数a,b,拟合多项式s(x)的系数c,其中d(k)=(y,sk); s1=0;s2=ones(1,m+1);v2=sum(w); d (1)=y*w';c (1)=d (1)/v2; fork=1: n xs=x.*s2.^2*w';a(k)=xs/v2; ifk==1b(k)=0; elseb(k)=v2/v1; end s3=(x-a(k)).*s2-b(k)*s1; v3=s3.^2*w'; d(k+1)=y.*s3*w';c(k+1)=d(k+1)/v3; s1=s2;s2=s3;v1=v2;v2=v3; end %求平方误差r r=y.*y*w'-c*d'; %求拟合多项式s(x)的降幂系数alph alph=zeros(1,n+1);T=zeros(n+1,n+2); T(: 2)=ones(n+1,1);T(2,3)=-a (1); ifn>=2 fork=3: n+1 fori=3: k+1 T(k,i)=T(k-1,i)-a(k-1)*T(k-1,i-1)-b(k-1)*T(k-2,i-2); end end end fori=1: n+1 fork=i: n+1 alph(n+2-i)=alph(n+2-i)+c(k)*T(k,k+2-i); end end %用秦九韶方法计算s(t)的输出序列(t,s) xmin=min(x);xmax=max(x);dx=(xmax-xmin)/(25*m); t=(xmin-dx): dx: (xmax+dx); s=alph (1); fork=2: n+1 s=s.*t+alph(k); end %输出点列x-y和拟合曲线t-s的图形 plot(x,y,'*',t,s,'-'); title('离散数据的多项式拟合'); xlabel('x');ylabel('y'); gridon; 实验3.1,输入1,权函数为[1111111],结果为: >>平方误差: 2.1762e-005 参数alph: 1.9991-2.9977-3.9683e-0050.54912 图形为: 实验3.2: 输入2,权函数为[1111111],结果为: >>平方误差: 2.1762e-005 参数alph: 1.9991-2.9977-3.9683e-0050.54912 图形为: 分析: 利用最小二乘法作曲线的拟合,对实验3.1给出的数据作的三次多项式的图形和数据节点拟合得较好。 正交化多项式,取w=[1111111]时,曲线的拟合也是十分吻合的。 由于没有出现病态法方程组的情况,二者图形结果及平方误差也相差不大。 第四章实验4 用matlab编写的程序如下: functiont_charpt4 %数值实验四: 含“实验4.1: 复化求积公式计算定积分”和“实验4.2: 高斯数值积分法用%于积分方程求解” %子函数调用: CG_L_Ifunction(复化Gauss_LegendreI)公式、CTrapezia(复化梯形%公式)、CSimpson(复化Simpson公式) %输入: 实验选择、积分法选择、积分式题号选择 %输出: 积分(实验4.1)或方程解(实验4.2)的精确值和数值解及误差 result=inputdlg({'请选择实验,若选4.1,请输入1,否则输入2: '},'charpt4',1,{'1'}); Nt=str2num(char(result)); if(Nt~=1)&(Nt~=2)errordlg('实验选择错误! ');return;end promps={'请选择积分法,若用复化梯形,输入T,用复化Simpson,输入S,用复化Gauss_Legendre,输入GL: '}; result=inputdlg(promps,'charpt4',1,{'T'}); Nb=char(result); if(Nb~='T'&Nb~='S'&Nb~='GL')errordlg('积分公式选择错误! ');return;end if(Nt==1) result=inputdlg({'请输入积分式题号1—4: '},'实验4.1',1,{'1'}); Nb_f=str2num(char(result)); if((Nb_f<1)|(Nb_f>4))errordlg('没有该积分式! ');return;end switchNb_f case1 fun=inline('-2./(x.^2-1)');a=2;b=3; case2 fun=inline('4./(x.^2+1)');a=0;b=1; case3 fun=inline('3.^x');a=0;b=1; case4 fun=inline('x.*exp(x)');a=1;b=2; end tol=0.5e-7;h=0.01; if(Nb=='T')%用复化梯形公式 t=(fun(a)+fun(b))*(b-a)/2; k=1;t0=0; while(abs(t-t0)>=tol*3) t0=t;h=(b-a)/2^k; t=t0/2+h*sum(fun(a+h: 2*h: b-h)); k=k+1; end elseif(Nb=='S')%用复化Simpson公式 t=quad(fun,a,b,tol); elseif(Nb=='GL')%用复化Gauss_LegendreI N=floor((b-a)/h);t=0;xk=0; fork=0: N xk=a+k*h+h/2; t=t+fun(xk-h/(2*sqrt(3)))+fun(xk+h/(2*sqrt(3))); end t=t*h/2; end elseif(Nt==2) result=inputdlg({'请输入方程式题号1或2: '},'实验4.2',1,{'1'}); Nb_f=str2num(char(result)); if(Nb_f~=1&Nb_f~=2)errordlg('没有该方程式! ');return;end result=inputdlg({'请输入步长: '},'实验4.2',1,{'0.01'}); h=str2num(char(result)); if(h<=0)errordlg('请输入正确的步长! ');return;end if(Nb=='T')%用复化梯形公式 [x,t]=CTrapezia(0,1,h,Nb_f); elseif(Nb=='S')%用复化Simpson公式 [x,t]=CSimpson(0,1,h,Nb_f); elseif(Nb=='GL')%用复化Gauss_LegendreI公式 [x,t]=CG_L_I(0,1,h,Nb_f); end plot(x,t,'g--'); xlabel('x');ylabel('y'); title('积分方程求解') holdon disp(['实验4.2(',num2str(Nb_f),')的计算结果: ',num2str(t')]); if(Nb_f==1) fplot('exp(x)',[01],'*'); holdoff disp(['实验4.2题 (1)的精确解: ',num2str(exp(x'))]); disp(['绝对误差和: ',num2str(sum(abs(t'-exp(x'))))]); else fplot('1/(1+t)^2',[01],'*'); holdoff y=1./(x+1).^2; disp(['实验4.2题 (2)的精确解: ',num2str(exp(y'))]); disp(['绝对误差和: ',num2str(sum(abs(t'-y')))]); end end if(Nt==1) disp(['实验4.1题(',num2str(Nb_f),')的计算结果: ',num2str(t)]); switchNb_f case1 disp('精确解: ln2-ln3=-0.4054651081') disp(['绝对误差: ',num2str(abs(t+0.4054651081))]); case2 disp('精确解: pi=3.14159265358979') disp(['绝对误差: ',num2str(abs(t-pi))]); case3 disp('精确解: 2/ln3=1.82047845325368') disp(['绝对误差: ',num2str(abs(t-1.82047845325368))]); case4 disp('精确解: e^2=7.38905609893065') disp(['绝对误差: ',num2str(abs(t-7.38905609893065))]); end end % function[x,y]=CG_L_I(a,b,h,N) %复化Gauss_LegendreI,用于积分方程求解 %输入: a、b分别为求积下、上限,h为步长,N为方程式题号实验选择、积分法选择、积分式题号选择 %输出: x,y为方程离散解y(i)=f(x(i)) n=floor((b-a)/h); A=zeros(2*n+2,2*n+2);b=zeros(2*n+2,1);x=b; fori=0: n t1=a+h*(i+0.5-sqrt(3)/6);t2=a+h*(i+0.5+sqrt(3)/6); x(2*i+1)=t1;x(2*i+2)=t2; if(N==1) b(2*i+1)=exp(t1);b(2*i+2)=exp(t2); else b(2*i+1)=-(4*t1.^3+5*t1.^2-2*t1+5)/(8*(t1+1).^2); b(2*i+2)=-(4*t2.^3+5*t2.^2-2*t2+5)/(8*(t2+1).^2); end forj=0: n x1=a+h*(j+0.5-sqrt(3)/6);x2=a+h*(j+0.5+sqrt(3)/6); if(N==1) A(2*i+1,2*j+1)=exp(t1)*2/(exp (1)-1); A(2*i+1,2*j+2)=exp(t1)*2/(exp (1)-1); A(2*i+2,2*j+1)=exp(t2)*2/(exp (1)-1); A(2*i+2,2*j+2)=exp(t2)*2/(exp (1)-1); else A(2*i+1,2*j+1)=1/(1+x1)-t1;A(2*i+1,2*j+2)=1/(1+x2)-t1; A(2*i+2,2*j+1)=1/(1+x1)-t2;A(2*i+2,2*j+2)=1/(1+x2)-t2; end end end A=h*A/2-eye(2*n+2);y=inv(A)*b; % function[x,y]=CTrapezia(a,b,h,N) %复化Gauss_%复化梯形公式,用于积分方程求解 %输入: a、b分别为求积下、上限,h为步长,N为方程式题号实验选择、积分法选择、积分式题号选择 %输出: x,y为方程离散解y(i)=f(x(i)) n=floor((b-a)/h); A=zeros(n+1,n+1);b=zeros(n+1,1);x=b; fori=0: n t=a+i*h;x(i+1)=t; if(N==1) b(i+1)=-exp(t); else b(i+1)=(4*t.^3+5*t.^2-2*t+5)/(8*(t+1).^2); end forj=0: n s=a+j*h; if(j==0)|(j==n) if(N==1) A(i+1,j+1)=exp(t)*2/(exp (1)-1); else A(i+1,j+1)=1/(1+s)-t; end else if(N==1) A(i+1,j+1)=exp(t)*2/(exp (1)-1)*2; else A(i+1,j+1)=(1/(1+s)-t)*2; end end end end A=A-eye(n+1)*2/h;b=-2*b/h;y=inv(A)*b; % function[x,y]=CSimpson(a,b,h,N) %复化Simpson公式,用于积分方程求解 %输入: a、b分别为求积下、上限,h为步长,N为方程式题号实验选择、积分法选择、积分式题号选择 %输出: x,y为方程离散解y(i)=f(x(i)) n=floor((b-a)/h);h1=h/2; A=zeros(2*n+1,2*n+1);b=zeros(2*n+1,1);x=b; fori=0: n*2 t=a+i*h1;x(i+1)=t; if(N==1) b(i+1)=-exp(t); else b(i+1)=(4*t.^3+5*t.^2-2*t+5)/(8*(t+1).^2); end forj=0: n*2 s=a+j*h1; if(j==0)|(j==n*2) if(N==1) A(i+1,j+1)=exp(t)*2/(exp (1)-1); else A(i+1,j+1)=1/(1+s)-t; end elseif(mod(j,2)==0) if(N==1) A(i+1,j+1)=exp(t)*2/(exp (1)-1)*2; else A(i+1,j+1)=(1/(1+s)-t)*2; end elseif(mod(j,2)~=0) if(N==1) A(i+1,j+1)=exp(t)*2/(exp (1)-1)*4; else A(i+1,j+1)=(1/(1+s)-t)*4; end end end end A=A-eye(2*n+1)*6/h;b=-6*b/h; y=inv(A)*b; 1.若选择实验4.1, 1)积分法选择T(即复化梯形公式): >>实验4.1题 (1)的计算结果: -0.40547 精确解: ln2-ln3=-0.4054651081 绝对误差: 1.3944e-008 >>实验4.1题 (2)的计算结果: 3.1416 精确解: pi=3.14159265358979 绝对误差: 3.9736e-008 >>实验4.1题(3)的计算结果: 1.8205 精确解: 2/ln3=1.82047845325368 绝对误差: 4.3655e-008 >>实验4.1题(4)的计算结果: 7.3891 精确解: e^2=7.38905609893065 绝对误差: 2.0775e-008 2)积分法选择S(即复化Simpson公式): >>实验4.1题 (1)的计算结果: -0.40547 精确解: ln2-ln3=-0.4054651081 绝对误差: 1.2625e-009 >>实验4.1题 (2)的计算结果: 3.1416 精确解: pi=3.14159265358979 绝对误差: 2.7517e-010 >>实验4.1题(3)的计算结果: 1.8205 精确解: 2/ln3=1.82047845325368 绝对误差: 1.0877e-010 >>实验4.1题(4)的计算结果: 7.3891 精确解: e^2=7.38905609893065 绝对误差: 8.2168e-011 3)积分法选择GL(即复化Gauss_Legendre公式): >>实验4.1题 (1)的计算结果: -0.40796 精确解: ln2-ln3=-0.4054651081 绝对误差: 0.0024907 >>实验4.1题 (2)的计算结果: 3.1615 精确解: pi=3
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 作业