数值分析实验报告Word文档下载推荐.docx
《数值分析实验报告Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《数值分析实验报告Word文档下载推荐.docx（17页珍藏版）》请在冰豆网上搜索。
end
result=inputdlg（{'
请输入插值结点数N:
},'
10'
Nd=str2num（char（result））;
if（Nd<
1）errordlg（'
结点输入错误！
switchNb_f
case'
f=inline（'
1./（1+25*x.^2）'
a=-1;
b=1;
x./（1+x.^4）'
a=-5;
b=5;
atan（x）'
b=5;
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:
if（j~=k）
p=p*（z-x0（j））/（x0（k）-x0（j））;
end
s=s+p*y0（k）;
y（i）=s;
实验结果分析
（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:
y=polyval（alph,x）;
plot（x,y,'
k--'
xlabel（'
y0*andpolyfit.y--'
holdon
plot（x0,y0,'
*'
）
gridon;
disp（['
平方误差:
num2str（r）]）
参数alph:
num2str（alph）]）
实验结果
2.1762e-005
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:
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%选取的点数
x=0:
1;
f=8.*（3*x.^2-1）./（x.^2+1）.^3;
x=1;
n=1/h
f=log（3）.*log（3）.*3.^x;
%plot（x,f）;
%画出二阶导函数图像
f=log（3）*log（3）*3^x;
x=1:
f=2.*exp（x）+x.*exp（x）;
%二阶导函数
x=2;
n=ceil（1/h）+1%选取的点数
估计结果步长h及结点数n分别为
h=0.8
n=1793
h=0.6
n=1827
h=0.9
n=2458
n=7020
（2）．用复化simpson公式进行事前估计的Matlab程序
f=-2*（（-72*x.^2-24）.*（x.^2-1）-192*x.^2.*（x.^2+1））./（x.^2-1）.^5;
%四阶导函数
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=2*ceil（1/h）+1%选取的点数
f=4*（（-72*x.^2+24）.*（x.^2+1）-192*x.^2.*（-x.^2+1））./（x.^2+1）.^5;
f=4*（（-72*x^2+24）*（x^2+1）-192*x^2*（-x^2+1））/（x^2+1）^5;
h=sqrt（sqrt（abs（h4）））%步长
f=log（3）^4*3.^x;
%计算导函数最大值
f=4*exp（x）+x.*exp（x）;
plot（x,f）%画出原函数
h=sqrt（sqrt（abs（h4）））
h=0.13411
n=47
h=0.76542
n=35
h=0.18433
n=29
h=0.18546
n=49
2.积分计算的Matlab程序：
promps={'
请选择积分公式，若用复化梯形，请输入T，用复化simpson，输入S，
用复化Gauss_Legendre，输入GL：
result=inputdlg（promps,'
charpt4'
T'
Nb=char（result）;
if（Nb~='
&
Nb~='
S'
GL'
）
errordlg（'
积分公式选择错误'
return;
end
result=inputdlg（{'
请输入积分式题号1-4：
实验4.1'
1'
Nb_f=str2num（char（result））;
if（Nb_f<
1|Nb_f>
4）
没有该积分式'
switchNb_f
case1
fun=inline（'
-2./（x.^2-1）'
a=2;
b=3;
case2
4./（x.^2+1）'
a=0;
b=1;
case3
3.^x'
case4
x.*exp（x）'
a=1;
b=2;
if（Nb=='
）%用复化梯形公式
promps={'
请输入用复化梯形公式应取的步长：
result=inputdlg（promps,'
实验4.2'
0.01'
h=str2num（char（result））;
if（h<
=0）
请输入正确的步长！
end
tic;
N=floor（（b-a）/h）;
detsum=0;
fori=1:
N-1
xk=a+i*h;
detsum=detsum+fun（xk）;
t=h*（fun（a）+fun（b）+2*detsum）/2;
time=toc;
t
）%用复化Simpson公式
请输入用复化Simpson公式应取的步长：
detsum_1=0;
detsum_2=0;
xk_1=a+i*h;
detsum_1=detsum_1+fun（xk_1）;
N
xk_2=a+h*（2*i-1）/2;
detsum_2=detsum_2+fun（xk_2）;
t=h*（fun（a）+fun（b）+2*detsum_1+4*detsum_2）/6;
）%用复化Gauss_LegendreI
%先根据复化Gauss_LegendreI公式的余项估计步长
请输入用复化Gauss_LegendreI公式应取的步长：
请输入正确的步长!
t=0;
fork=0:
xk=a+k*h+h/2;
t=t+fun（xk-h/（2*sqrt（3）））+fun（xk+h/（2*sqrt（3）））;
t=t*h/2;
disp（'
精确解：
ln2-ln3=-0.4054651081'
disp（['
绝对误差：
num2str（abs（t+0.4054651081））]）;
运行时间：
num2str（time）]）;
pi=3.979'
num2str（abs（t-pi））]）;
2/ln3=1.368'
num2str（abs（t-1.368））]）;
e^2=7.065'
num2str（abs（t-7.065））]）;
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公式进行计算时：
t=-0.7519
2.7519e-011
0.022
t=3.979
0
0.021
t=1.288
9.2018e-012
0.019
t=7.118
9.0528e-011
3.当选用复化Gauss-LegendreI型公式进行计算时：
t=-0.5262
4.7385e-012
0.023
1.3323e-015
t=1.754
6.1431e-012
t=7.046
1.441e-012
结果分析：
当选用复化梯形公式时，对步长的事前估计所要求的步长很小，选取的节点很多，误差绝对限要达到
时，对不同的函数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;
请输入[-50，50]间的参数a:
实验5.1'
-40'
a=str2num（char（result））;
if（a<
-50|a>
50）errordlg（'
请输入正确的参数a!
请输入（01）之间的步长:
h=str2num（char（result））;
if（h<
=0|h>
=1）errordlg（'
请输入正确的（01）间的步长!
h:
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;
升级会员 