数值分析实验报告.docx
《数值分析实验报告.docx》由会员分享，可在线阅读，更多相关《数值分析实验报告.docx（17页珍藏版）》请在冰豆网上搜索。
数值分析实验报告
实验2.1多项式插值的振荡现象
实验目的：
在一个固定的区间上用插值逼近一个函数，显然Lagrange插值中使用的节点越多，插值多项式的次数就越高。
我们自然关心插值多项式的次数增加时，Ln（x）是否也更加靠近被逼近的函数。
Runge给出的一个例子是极著名并富有启发性的。
实验容：
设区间[-1，1]上函数f（x）=1/（1+25x2）。
考虑区间[-1,1]的一个等距划分，分点为xi=-1+2i/n，i=0，1，2，…，n，
则拉格朗日插值多项式为
.
其中，li（x），i=0，1，2，…，n是n次Lagrange插值基函数。
实验步骤与结果分析：
实验源程序
functionChap2Interpolation
%数值实验二：
“实验2.1：
多项式插值的震荡现象”
%输入：
函数式选择，插值结点数
%输出：
拟合函数及原函数的图形
promps={'请选择实验函数，若选f（x）,请输入f,若选h（x）,请输入h,若选g（x）,请输入g:
'};
titles='charpt_2';
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:
'},'charpt_2',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
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）--'）;
%--------------------------------------------------------------------
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
实验结果分析
（1）增大分点n=2，3，…时，拉格朗日插值函数曲线如图所示。
n=6
n=7
n=8
n=9
n=10
从图中可以看出，随着n的增大，拉格朗日插值函数在x=0附近较好地逼近了原来的函数f（x），但是却在两端x=-1和x=1处出现了很大的振荡现象。
并且，仔细分析图形，可以看出，当n为奇数时，虽然有振荡，但振荡的幅度不算太大，n为偶数时，其振荡幅度变得很大。
通过思考分析，我认为，可能的原因是f（x）本身是偶函数，如果n为奇数，那么Lagrange插值函数Ln（x）的最高次项xn-1是偶次幂，比较符合f（x）本身是偶函数的性质；如果n为偶数，那么Lagrange插值函数Ln（x）的最高次项xn-1是奇次幂，与f（x）本身是偶函数的性质相反，因此振荡可能更剧烈。
（2）将原来的f（x）换为其他函数如h（x）、g（x），结果如图所示。
其中h（x）,g（x）均定义在[-5，5]区间上，h（x）=x/（1+x4），g（x）=arctanx。
h（x）,n=7
h（x）,n=8
h（x）,n=9
h（x）,n=10
g（x）,n=7
g（x）,n=8
g（x）,n=9
g（x）,n=10
分析两个函数的插值图形，可以看出：
随着n的增大，拉格朗日插值函数在x=0附近较好地逼近了原来的函数f（x），但是却在两端x=-5和x=5处出现了很大的振荡现象。
并且，仔细分析图形，可以看出，当n为偶数时，虽然有振荡，但振荡的幅度不算太大，n为奇数时，其振荡幅度变得很大。
原因和上面f（x）的插值类似，h（x）、g（x）本身是奇函数，如果n为偶数，那么Lagrange插值函数Ln（x）的最高次项xn-1是奇次幂，比较符合h（x）、g（x）本身是奇函数的性质；如果n为奇数，那么Lagrange插值函数Ln（x）的最高次项xn-1是偶次幂，与h（x）、g（x）本身是奇函数的性质相反，因此振荡可能更剧烈。
实验3.1多项式最小二乘拟合
实验目的：
编制以函数{xk}k=0,…,n；为基的多项式最小二乘拟合程序。
实验容：
对表中的数据作三次多项式最小二乘拟合。
xi
-1.0
-0.5
0.0
0.5
1.0
1.5
2.0
yi
-4.447
-0.452
0.551
0.048
-0.447
0.549
4.552
取权函数wi≡1，求拟合曲线
中的参数{αk}、平方误差δ2，并作离散据
{xi,yi}的拟合函数的图形。
实验源程序
functionChap3CurveFitting
%数值实验三:
“实验3.1”
%输出:
原函数及求得的相应插值多项式的函数的图像以及参数alph和误差r
x0=-1:
0.5:
2;
y0=[-4.447-0.4520.5510.048-0.4470.5494.552];
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*andpolyfit.y--'）;
holdon
plot（x0,y0,'*'）
gridon;
disp（['平方误差:
',num2str（r）]）
disp（['参数alph:
',num2str（alph）]）
实验结果
平方误差:
2.1762e-005
参数alph:
1.9991-2.9977-3.9683e-0050.54912
实验4.1
实验目的：
复化求积公式计算定积分.
实验题目：
数值计算下列各式右端定积分的近似值.
实验要求：
（1）若用复化梯形公式、复化Simpson公式和复化Gauss-LegendreI型公式做计算，要求绝对误差限为
，分别利用它们的余项对每种算法做出步长的事前估计.
（2）分别用复化梯形公式，复化Simpson公式和复化Gauss-LegendreI型公式作计算.
（3）将计算结果与精确解做比较，并比较各种算法的计算量.
实验程序：
1.事前估计的Matlab程序如下：
（1）．用复化梯形公式进行事前估计的Matlab程序
formatlongg
x=2:
0.01:
3;
f=-4*（3*x.^2+1）./（x.^2-1）.^3;%二阶导函数
%plot（x,f）%画出二阶导函数图像
x=2.0;%计算导函数最大值
f=-4*（3*x^2+1）/（x^2-1）^3;
h2=0.5*10^（-7）*12/f;
h=sqrt（abs（h2））%步长
n=1/h;
n=ceil（1/h）+1%选取的点数
formatlongg
x=0:
0.01:
1;
f=8.*（3*x.^2-1）./（x.^2+1）.^3;%二阶导函数
%plot（x,f）%画出二阶导函数图像
x=1;%计算导函数最大值
f=8.*（3*x.^2-1）./（x.^2+1）.^3;
h2=0.5*10^（-7）*12/f;
h=sqrt（abs（h2））%步长
n=1/h
n=ceil（1/h）+1%选取的点数
formatlongg
x=0:
0.01:
1;
f=log（3）.*log（3）.*3.^x;%二阶导函数
%plot（x,f）;%画出二阶导函数图像
x=1;%计算导函数最大值
f=log（3）*log（3）*3^x;
h2=0.5*10^（-7）*12/f;
h=sqrt（abs（h2））%步长
n=1/h
n=ceil（1/h）+1%选取的点数
formatlongg
x=1:
0.01:
2;
f=2.*exp（x）+x.*exp（x）;%二阶导函数
%plot（x,f）%画出二阶导函数图像
x=2;%计算导函数最大值
f=2.*exp（x）+x.*exp（x）;
h2=0.5*10^（-7）*12/f;
h=sqrt（abs（h2））%步长
n=1/h
n=ceil（1/h）+1%选取的点数
估计结果步长h及结点数n分别为
h=0.8
n=1793
h=0.6
n=1827
h=0.9
n=2458
h=0.9
n=7020
（2）．用复化simpson公式进行事前估计的Matlab程序
formatlongg
x=2:
0.01:
3;
f=-2*（（-72*x.^2-24）.*（x.^2-1）-192*x.^2.*（x.^2+1））./（x.^2-1）.^5;%四阶导函数
x=2.0;
f=-2*（（-72*x^2-24）*（x^2-1）-192*x^2*（x^2+1））/（x^2-1）^5;%计算导函数最大值
h4=0.5*10^（-7）*180*16/f;
h=sqrt（sqrt（abs（h4）））%步长
n=1/h;%求分段区间个数
n=2*ceil（1/h）+1%选取的点数
formatlongg
x=0:
0.01:
1;
f=4*（（-72*x.^2+24）.*（x.^2+1）-192*x.^2.*（-x.^2+1））./（x.^2+1）.^5;%四阶导函数
x=1;
f=4*（（-72*x^2+24）*（x^2+1）-192*x^2*（-x^2+1））/（x^2+1）^5;%计算导函数最大值
h4=0.5*10^（-7）*180*16/f;
h=sqrt（sqrt（abs（h4）））%步长
n=1/h;%求分段区间个数
n=2*ceil（1/h）+1%选取的点数
formatlongg
x=0:
0.01:
1;
f=log（3）^4*3.^x;%四阶导函数
x=1;
f=log（3）^4*3.^x;%计算导函数最大值
h4=0.5*10^（-7）*180*16/f;
h=sqrt（sqrt（abs（h4）））%步长
n=1/h;%求分段区间个数
n=2*ceil（1/h）+1%选取的点数
formatlongg
x=1:
0.01:
2;
f=4*exp（x）+x.*exp（x）;%四阶导函数
plot（x,f）%画出原函数
x=2;
f=4*exp（x）+x.*exp（x）;%计算导函数最大值
h4=0.5*10^（-7）*180*16/f;
h=sqrt（sqrt（abs（h4）））
n=1/h;%求分段区间个数
n=2*ceil（1/h）+1%选取的点数
估计结果步长h及结点数n分别为
h=0.13411
n=47
h=0.76542
n=35
h=0.18433
n=29
h=0.18546
n=49
2.积分计算的Matlab程序：
formatlongg
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
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
if（Nb=='T'）%用复化梯形公式
promps={'请输入用复化梯形公式应取的步长：
'};
result=inputdlg（promps,'实验4.2',1,{'0.01'}）;
h=str2num（char（result））;
if（h<=0）
errordlg（'请输入正确的步长！
'）;
return;
end
tic;
N=floor（（b-a）/h）;
detsum=0;
fori=1:
N-1
xk=a+i*h;
detsum=detsum+fun（xk）;
end
t=h*（fun（a）+fun（b）+2*detsum）/2;
time=toc;
t
end
if（Nb=='S'）%用复化Simpson公式
promps={'请输入用复化Simpson公式应取的步长：
'};
result=inputdlg（promps,'实验4.2',1,{'0.01'}）;
h=str2num（char（result））;
if（h<=0）
errordlg（'请输入正确的步长！
'）;
return;
end
tic;
N=floor（（b-a）/h）;
detsum_1=0;
detsum_2=0;
fori=1:
N-1
xk_1=a+i*h;
detsum_1=detsum_1+fun（xk_1）;
end
fori=1:
N
xk_2=a+h*（2*i-1）/2;
detsum_2=detsum_2+fun（xk_2）;
end
t=h*（fun（a）+fun（b）+2*detsum_1+4*detsum_2）/6;
time=toc;
t
end
if（Nb=='GL'）%用复化Gauss_LegendreI
%先根据复化Gauss_LegendreI公式的余项估计步长
promps={'请输入用复化Gauss_LegendreI公式应取的步长：
'};
result=inputdlg（promps,'实验4.2',1,{'0.01'}）;
h=str2num（char（result））;
if（h<=0）
errordlg（'请输入正确的步长!
'）;
return;
end
tic;
N=floor（（b-a）/h）;t=0;
fork=0:
N-1
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;
time=toc;
t
end
switchNb_f
case1
disp（'精确解：
ln2-ln3=-0.4054651081'）
disp（['绝对误差：
',num2str（abs（t+0.4054651081））]）;
disp（['运行时间：
',num2str（time）]）;
case2
disp（'精确解：
pi=3.979'）
disp（['绝对误差：
',num2str（abs（t-pi））]）;
disp（['运行时间：
',num2str（time）]）;
case3
disp（'精确解：
2/ln3=1.368'）
disp（['绝对误差：
',num2str（abs（t-1.368））]）;
disp（['运行时间：
',num2str（time）]）;
case4
disp（'精确解：
e^2=7.065'）
disp（['绝对误差：
',num2str（abs（t-7.065））]）;
disp（['运行时间：
',num2str（time）]）;
end
1.当选用复化梯形公式时：
（1）式运行结果为：
t=-0.351
精确解：
ln2-ln3=-0.4054651081
绝对误差：
1.3944e-008
运行时间：
0.003
（2）式运行结果为：
t=3.336
精确解：
pi=3.979
绝对误差：
3.9736e-008
运行时间：
0.005
（3）式运行结果为：
t=1.861
精确解：
2/ln3=1.368
绝对误差：
4.3655e-008
运行时间：
0.016
（4）式运行结果为：
t=7.610
精确解：
e^2=7.065
绝对误差：
2.0775e-008
运行时间：
0.007
2.当选用复化Simpson公式进行计算时：
（1）式运行结果为：
t=-0.7519
精确解：
ln2-ln3=-0.4054651081
绝对误差：
2.7519e-011
运行时间：
0.022
（2）式运行结果为：
t=3.979
精确解：
pi=3.979
绝对误差：
0
运行时间：
0.021
（3）式运行结果为：
t=1.288
精确解：
2/ln3=1.368
绝对误差：
9.2018e-012
运行时间：
0.019
（4）式运行结果为：
t=7.118
精确解：
e^2=7.065
绝对误差：
9.0528e-011
运行时间：
0.021
3.当选用复化Gauss-LegendreI型公式进行计算时：
（1）式运行结果为：
t=-0.5262
精确解：
ln2-ln3=-0.4054651081
绝对误差：
4.7385e-012
运行时间：
0.023
（2）式运行结果为：
t=3.979
精确解：
pi=3.979
绝对误差：
1.3323e-015
运行时间：
0.021
（3）式运行结果为：
t=1.754
精确解：
2/ln3=1.368
绝对误差：
6.1431e-012
运行时间：
0.019
（4）式运行结果为：
t=7.046
精确解：
e^2=7.065
绝对误差：
1.441e-012
运行时间：
0.021
结果分析：
当选用复化梯形公式时，对步长的事前估计所要求的步长很小，选取的节点很多，误差绝对限要达到
时，对不同的函数n的取值需达到1000-10000之间，计算量是很大。
用复化simpson公式对步长的事前估计所要求的步长相对大些，选取的节点较少，误差绝对限要达到
时，对不同的函数n的取值只需在10-100之间，计算量相对小了很多，可满足用较少的节点达到较高的精度，比复化梯形公式的计算量小了很多。
用复化simpson公式计算所得的结果比用复化梯形公式计算所得的结果精度高很多，而且计算量小。
当选用Gauss-LagrangeI型公式进行计算时，选用较少的节点就可以达到很高的精度。
实验5.1常微分方程性态和R-K法稳定性试验
实验目的：
考察下面微分方程右端项中函数y前面的参数对方程性态的影响（它可使方程为好条件的或坏条件的）和研究计算步长对R-K法计算稳定性的影响。
实验容及要求：
实验题目：
常微分方程初值问题
其中，
。
其精确解为
。
实验要求：
本实验题都用4阶经典R-K法计算。
（1）对参数a分别取4个不同的数值：
一个大的正值，一个小的正值，一个绝对值小的负值和一个绝对值大的负值。
取步长h=0.01，分别用经典的R-K法计算，将四组计算结果画在同一图上，进行比较并说明相应初值问题的性态。
（2）取参数a为一个绝对值不大的负值和两个计算步长，一个步长使参数ah在经典R-K法的稳定域，另一个步长在经典R-K法的稳定域外。
分别用经典R-K法计算并比较计算结果。
取全域等距的10个点上的计算值，列表说明。
实验程序：
Matlab程序如下：
functioncharp5RK
%数值试验5.1：
常微分方程性态和R-K法稳定性试验
%输入：
参数a，步长h
%输出：
精确解和数值解图形对比
%clf;
result=inputdlg（{'请输入[-50，50]间的参数a:
'},'实验5.1',1,{'-40'}）;
a=str2num（char（result））;
if（a<-50|a>50）errordlg（'请输入正确的参数a!
'）;return;end
result=inputdlg（{'请输入（01）之间的步长:
'},'实验5.1',1,{'0.01'}）;
h=str2num（char（result））;
if（h<=0|h>=1）errordlg（'请输入正确的（01）间的步长!
'）;return;end
x=0:
h:
1;
y=x;
N=length（x）;
y
（1）=1;
func=inline（'1+（y-x）.*a'）;
forn=1:
N-1
k1=func（a,x（n）,y（n））;
k2=func（a,x（n）+h/2,y（n）+k1*h/2）;
k3=func（a,x（n）+h/2,y（n）+k2*h/2）;
k4=func（a,x（n）+h,y（n）+k3*h）;
y（n+1）=y（n）+h*（k1+2*k2+2*k3+k4）/6;
升级会员 