数值分析实验题.docx
- 文档编号:17112457
- 上传时间:2023-07-22
- 格式:DOCX
- 页数:27
- 大小:572.84KB
数值分析实验题.docx
《数值分析实验题.docx》由会员分享,可在线阅读,更多相关《数值分析实验题.docx(27页珍藏版)》请在冰点文库上搜索。
数值分析实验题
数值分析实验报告
第一题实验题1.2
1、实验内容
实验1.2体会稳定性在选择算法中的地位,误差扩张的算法不稳定,而误差衰竭的算法是稳定的。
分别采用E.1.6(即E.1.4)
和算法算法E.1.7(即E.1.5)
两种算法。
2、源程序
%functiont_charpt1
%数值试验1.2:
误差传播与算法稳定性
%输入:
递推公式选择与递推步数
%输出:
各步递推值及误差结果,以及递推值和误差与递推步数的关系图
clear;
clc;
promps={'请选择递推关系式,若选(1.4),请输人1,否则输入2:
'};
result=inputdlg(promps,'charpt1_2',1,{'1'});
Nb=str2num(char(result));
if((Nb~=1)&(Nb~=2))
errordlg('请选择递推关系式,若选(1.4),请输人1,否则输人2!
');
return;
end
result=inputdlg({'请输人递推步数n:
'},'charpt1_2',1,{'10'});
steps=str2num(char(result));
if(steps<1)
errordlg('递推步数错误!
');
return;
end
result=inputdlg({'请输入计算中所采用的有效数字位数:
'},'charpt1_2',1,{'5'});
Sd=str2num(char(result));
formatlong%设置显示精度
result=zeros(1,steps);%存储计算结果
err=result;%存储计算的绝对误差
func=result;%存储用库函数quadl计算出的积分的近似值
%用库函数quadl计算积分的近似值
forn=1:
steps
fun=@(x)x.^n.*exp(x-1);
func(n)=quadl(fun,0,1);
end
if(Nb==1)
%用算法(1.4)计算
digits(Sd);%控制有效数字位数
result
(1)=subs(vpa(1/exp
(1)));
forn=2:
1:
steps
result(n)=subs(vpa(1-n*result(n-1)));
end
err=abs(result-func);
elseif(Nb==2)
%用算法(1.5)计算
digits(Sd);%控制有效数字位数
result(steps)=0;
forn=steps:
-1:
2
result(n-1)=subs(vpa((1-result(n))/n));
end
err=abs(result-func);
end
clf;%清除当前图像窗口
disp('递推值:
');
disp(sprintf('%e',result));
disp('误差:
');
disp(sprintf('%e',err));
plot([1:
steps],result,'-','LineWidth',2);
set(gca,'linewidth',0.5,'fontsize',16);
gridon
holdon;
plot([1:
steps],err,'r--','LineWidth',2);
xlabel('stepsn','FontSize',18);
ylabel('En-andERRn--','FontSize',18);
legend('En','err(n)');
title(['Algorithm(1.',num2str(Nb+3),')SignificantDigits',num2str(Sd)],'FontSize',18);
%text(2,err
(2),'\uparrowerr(n)');
%text(4,result(4),'\downarrowEn');
3、实验结果
(1)算法E1.6,有效数字5位
递推值:
3.678800e-012.642400e-012.072800e-011.708800e-011.456000e-011.264000e-011.152000e-017.840000e-022.944000e-01-1.944000e+00
误差:
5.588280e-071.117662e-063.352927e-061.341222e-056.705713e-054.023702e-042.816427e-032.253226e-022.027877e-012.02
(2)算法E1.6,有效数字6位
递推值:
3.678790e-012.642420e-012.072740e-011.709040e-011.454800e-011.271200e-011.101600e-011.187200e-01-6.848000e-021.684800e+00
误差:
4.411720e-078.823378e-072.647073e-061.058778e-055.294287e-053.176298e-042.223573e-031.778774e-021.600923e-011.60
(3)算法E1.6,有效数字7位
递推值:
3.678794e-012.642412e-012.072764e-011.708944e-011.455280e-011.268320e-011.121760e-011.025920e-017.667200e-022.332800e-01
误差:
4.117197e-088.233779e-082.470726e-079.877761e-074.942873e-062.962984e-052.075730e-041.659738e-031.494029e-021.49
(4)算法E1.7,有效数字5位
递推值:
3.678800e-012.642400e-012.072800e-011.708900e-011.455300e-011.267900e-011.125000e-011.000000e-011.000000e-010.000000e+00
误差:
5.588280e-071.117662e-063.352927e-063.412224e-062.942873e-061.237016e-051.164270e-049.322618e-048.387707e-038.38
(5)算法E1.7,有效数字6位
递推值:
3.678800e-012.642410e-012.072770e-011.708930e-011.455360e-011.267860e-011.125000e-011.000000e-011.000000e-010.000000e+00
误差:
5.588280e-071.176622e-073.529274e-074.122239e-073.057127e-061.637016e-051.164270e-049.322618e-048.387707e-038.38
(6)算法E1.7,有效数字7位
递推值:
3.678795e-012.642411e-012.072768e-011.708929e-011.455357e-011.267857e-011.125000e-011.000000e-011.000000e-010.000000e+00
误差:
5.882803e-081.766221e-081.529274e-075.122239e-072.757127e-061.667016e-051.164270e-049.322618e-048.387707e-038.38
4、结果分析
采用算法E1.7(即算法E1.5)能得到更精确的结果,当然,有效数字越多,结果越准确。
当采用算法E.1.6(即算法E.1.4)时:
设
的真实值为
则真实值
(1)
又有
(2)
(1)
(2)得:
(3)
对(3)式两边取绝对值得
(4)
由(4)可计算得
(5)
当采用算法E.1.7(即算法E.1.5)时:
同理:
设
的真实值为
则真实值
(6)
又有
(7)
(6)
(7)得:
(8)
对(8)式两边取绝对值得
(9)
由(9)可计算得
(10)
算法E.1.6(即算法E.1.4)的
很小,当n增大时,
增长速度很快,而算法E.1.7(即算法E.1.5)中的
虽然比较大,但是当n减小时,
呈现递减趋势。
所以比较两个算法,当某一步产生误差后,算法E.1.6(即算法E.1.4)对后面的影响是扩张的,而算法E.1.7(即算法E.1.5)对后面的影响是衰减的。
通过理论分析与计算实验,算法E1.7(即算法E1.5)更加稳定。
第二题实验题3.1
1、实验内容
实验3.1编制以函数
为基的多项式最小二乘拟合程序
2、源程序
functioncharpt3
formatlong;
%数值实验三:
含"实验3.1"和"实验3.2"
%子函数调用:
dlsa
%输入:
实验选择
%输出:
原函数及求得的相应插值多项式的函数的图像以及参数alph和误差r
result=inputdlg({'请选择实验,若选3.1,请输入1,否则输入2:
'},'charpt_3',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];
n=3;%n为拟合阶次
if(Nb==1)
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('y');
holdon;
plot(x0,y0,'*');
title('离散数据的多项式拟合');
gridon;
else
result=inputdlg({'请输入权向量w:
'},'charpt_3',1,{'[1111111]'});
w=str2num(char(result));
[a,b,c,alph,r]=dlsa(x0,y0,w,n);
end
disp(['平方误差:
',sprintf('%g',r)]);
disp(['参数alph:
',sprintf('%g\t',alph)])
%-------------------------------------------------------------------------
function[a,b,c,alph,r]=dlsa(x,y,w,n)
%功能:
用正交化方法对离散数据作多项式最小二乘拟合。
%输入:
m+1个离散点(x,y,w),x,y,w分别用行向量给出。
%拟合多项式的次数n,0 %输出: 三项递推公式的参数a,b,拟合多项式s(x)的系数c和alph, %平方误差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; if(k==1) b(k)=0; else b(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); if(n>=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、实验结果 平方误差: 2.17619e-05 参数alph: 1.99911-2.99767-3.96825e-050.549119 4、结果分析 利用最小二乘法作曲线的拟合,对实验3.1给出的数据作的三次多项式的图形和数据节点拟合得较好。 平方误差的数量级是10的-5次方,因此拟合效果很好。 最小二乘曲线拟合实际上是在离散情形下的最佳平方逼近.由于实验观测数据本身带有误差,所求的近似曲线并不要求过所有的给定点,只要求逼近函数能反映数据的变化趋势,最好的标准是要求观测函数 与 的偏差[ - ]的平方和 最小。 第三题实验题4.1 1、实验内容 实验4.1复化求积公式计算定积分 2、源程序 functioncharpt4 %数值实验四: 含“实验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为步
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 实验