数值分析实验题.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为步
升级会员 