南京邮电大学 数值计算实践报告.docx
《南京邮电大学 数值计算实践报告.docx》由会员分享，可在线阅读，更多相关《南京邮电大学 数值计算实践报告.docx（46页珍藏版）》请在冰点文库上搜索。
南京邮电大学数值计算实践报告
数值计算实践
、方程求根
1、实验目的
熟悉和掌握Newton法,割线法,抛物线法的方法思路，并能够在matlab上编程实现
2、问题描述
（1）.给定一个三次方程,分别用Newton法,割线法,抛物线法求解.
方程的构造方法:
（a）根:
方程的根为学号的后三位乘以倒数第二位加1再除以1000.
假设你的学号为B06060141,则根为141*（4+1）/1000=0.564
（b）方程:
以你的学号的后三位数分别作为方程的三次项,二次项,一次项的系数,根据所给的根以及三个系数确定常数项.
例如:
你的学号是B06060141,则你的方程是x3+4x2+x+a0=0的形式.
方程的根为0.564,因此有
0.5643+4*0.5642+0.564+a0=0,于是a0=-2.015790144
你的方程为x3+4x2+x-2.015790144=0.
（2）假设方程是sinx+4x2+x+a0=0的形式（三个系数分别是学号中的数字），重新解决类似的问题
（3）构造一个五次方程完成上面的工作.
四次方程的构造:
将三次多项式再乘以（x-p*）2得到对应的五次多项式（p*为已经确定的方程的根,显然,得到的五次方程有重根）.
（4）将
（2）中的方程同样乘以（x-p*）得到一个新的方程来求解
注:
（1）Newton法取0.5为初值,割线法以0,1为初值,抛物线法以0,0.5,1为初值,
（2）计算精度尽量地取高.
终止准则：
根据
来终止
（3）可供研究的问题：
（一）
的取值不同对收敛速度有多大的影响
（二）将注
（1）中的初值该为其它的初值，对收敛性以及收敛速度有无影响
（三）能否求出方程的所有的根
（4）实验报告的撰写
实验报告包含的内容：
（一）实验目的
（二）问题描述（三）算法介绍（包括基本原理）（四）程序（五）计算结果（六）结果分析（七）心得体会
3、算法介绍
在本问题中，我们用到了newton法，割线法，抛物线法。
1.Newton法迭代格式为：
当初值与真解足够靠近，newton迭代法收敛，对于单根，newton收敛速度很快，对于重根，收敛较慢。
2.割线法：
为了回避导数值的计算，使用上的差商代替，得到割线法迭代公式：
割线法的收敛阶虽然低于newton法，但迭代以此只需计算一次函数值，不需计算其导数，所以效率高，实际问题中经常应用。
3.抛物线法：
可以通过三点做一条抛物线，产生迭代序列的方法称为抛物线法。
其迭代公式为：
其中
是一阶均差和二阶均差。
收敛速度比割线法更接近于newton法。
对于本问题的解决就以上述理论为依据。
终止准则为：
本题中所有精度取1e-8。
4、程序计算结果
问题一
根据所给的要求，可知待求的方程为：
牛顿法
建立newton_1.m的源程序，源程序代码为：
functiony=newton_1（a,n,x0,nn,eps1）
x
（1）=x0;
b=1;
i=1;
while（abs（b）>eps1*x（i））
i=i+1;
x（i）=x（i-1）-n_f（a,n,x（i-1））/n_df（a,n,x（i-1））;
b=x（i）-x（i-1）;
if（i>nn）error
return;
end
end
y=x（i）;
建立n_f.m的源程序来求待求根的实数代数方程的函数，源程序代码为：
functiony=n_f（a,n,x）%待求根的实数代数方程的函数
y=0.0;
fori=1:
（n+1）
y=y+a（i）*x^（n+1-i）;
end
建立n_df.m的源程序来方程的一阶导数的函数，源程序代码为：
functiony=n_df（a,n,x）%方程的一阶导数的函数
y=0.0;
fori=1:
n
y=y+a（i）*（n+1-i）*x^（n-i）;
end
在matlab软件中执行下列语句并得到的最终结果截图：
割线法
建立gexian.m的源程序，源程序代码为
functionx=gexian（f,x0,x1,e）
ifnargin<4,e=1e-4;end
y=x0;x=x1;i=0;
whileabs（x-y）>e
i=i+1;
z=x-（feval（f,x）*（x-y））/（feval（f,x）-feval（f,y））;
y=x;
x=z;
end
i
在matlab软件中执行下列语句并得到的最终结果截图：
抛物线
建立paowuxian.m的源程序，源程序代码为：
functionx=paowuxian（f,x0,x1,x2,e）
ifnargin<4,e=1e-4;end
x=x2;y=x1;z=x0;i=0;
whileabs（x-y）>e
i=i+1;
h1=y-z;
h2=x-y;
c1=（feval（f,y）-feval（f,z））/h1;
c2=（feval（f,x）-feval（f,y））/h2;
d=（c1-c2）/（h1+h2）;
w=c2+h2*d;
xi=x-（2*feval（f,x））/（w+（w/abs（w））*sqrt（w^2-4*feval（f,x）*d））;
z=y;
y=x;
x=xi;
end
i
在matlab软件中执行下列语句并得到的最终结果截图：
研究一：
只改变初值
由上述结果可知，方程的解在0.2附近，所以将牛顿法为0.2；割线法的初值设为0，0.4；抛物线法的初值设为0，0.2，0.4；
牛顿法
根据问题1中牛顿法的程序，在matlab软件中执行下列语句并得到的最终结果截图：
割线法
根据问题1中割线程序，在matlab软件中执行下列语句并得到的最终结果截图：
抛物线法
根据问题1中抛物线法程序，在matlab软件中执行下列语句并得到的最终结果截图：
研究二只改变精度
将精度由1e-8改为1e-50和1e-100观察迭代次数有何变化
牛顿法：
根据问题1中的牛顿法的程序，在matlab软件中执行下列语句并得到的最终结果截图：
精度为1e-50时
精度为1e-100时
割线法
根据问题1中的割线法的程序，在matlab软件中执行下列语句并得到的最终结果截图：
精度为1e-50时
精度为1e-100时
抛物线法
根据问题1中的抛物线法的程序，在matlab软件中执行下列语句并得到的最终结果截图：
精度为1e-50时
精度为1e-100时
、
研究结论
在只改变初值时，当初值定得越靠近初值，迭代次数就越少。
在只改变精度时，当精度越来越大时，迭代次数并几乎不变。
综上所述，初值对迭代次数的影响比较大，精度对迭代次数影响不大。
问题二
问题描述
根据所给的要求，可知待求的方程为：
问题分析
仍然利用
（1）中方法求解这一问题，并利用图解法找到初值，通过观察图像，将newton法初值设为：
0.1，割线法初值设为：
0,0.2。
抛物线法初值设为：
0,0.1,0.2。
图像见下图：
Newton法
问题一的牛顿法的求解只适用于线性方程，所以在问题二中用其他方法来求解方程。
建立newton1.m源程序，源程序代码为：
functionx=newton1（fn,dfn,x0,e）
ifnargin<4,e=1e-4;end
x=x0;x0=x+2*e;i=1；
whileabs（x0-x）>e
i=i+1；
x0=x;
x=x0-feval（fn,x0）/feval（dfn,x0）;
end
在matlab软件中执行下列语句并得到最终结果截图
割线法
利用问题一中的割线法程序，在matlab软件中执行下列语句
抛物线法
利用问题一中的抛物线法程序，在matlab软件中执行下列语句
问题三
问题描述
按照题目要求对五次方程进行构造为：
问题分析
仍然利用一中方法求解这一问题，并利用图解法找到初值，通过观察图像，将牛顿法的两组初值为0以及0.5，割线法初值设为：
0,0.5以及0，0.2。
抛物线法初值设为：
两组初值为0,0.5,1以及0,0.25,0.5。
牛顿法
利用问题二中的程序，在matlab软件中执行下列语句：
初值为0时
初值为0.5时
割线法
利用问题一中割线法的程序，在matlab软件中执行下列语句：
初值为0，0.5时
初值为0，0.2时
抛物线法
利用问题一中抛物线法的程序，在matlab软件中执行下列语句：
初值为0，0.5，1时
初值为0，0.25，0.5时
问题四
问题描述
根据题目要求对方程进行构造为：
问题分析
仍然利用问题一中方法求解这一问题，并利用图解法找到初值，通过观察图像，newton法初值选取了两组初值为0以及0.5，割线法初值设为：
0，0.5和0，0.3。
抛物线法初值设为：
两组初值为0,0.2,0.4以及0，0.5，1。
牛顿法
利用问题二中的程序，在matlab软件中执行下列语句：
初值为0时
初值为0.5时
割线法
利用问题一中割线法的程序，在matlab软件中执行下列语句：
初值为0，0.5时
初值为0，0.3时
抛物线法
利用问题一中抛物线法的程序，在matlab软件中执行下列语句：
初值为0，0.1，0.2时
初值为0，0.5，1时
5、计算结果及分析
Newton法
割线法
抛物线法
问题一
0.224
0.224
0.224
问题二
0.1466
0.1466
0.1466
问题三
0.224
0.224
0.224
问题四
0.224和
0.146577258557101
0.224和
0.1466
0.224和
0.1466
问题一
迭代步数
问题二
迭代步数
问题三
迭代步数
问题四
迭代步数
Newton法
6
5
37
89
割线法
8
5
50
1310
抛物线法
9
5
63
1413
结果分析
将Newton法，割线法，抛物线法进行比较可以看到在本文题中，三种方法计算得到的最终结果基本相同，但是迭代步数有较大差别，综合看来Newton法迭代步数最少，割线法次之，抛物线法最次。
在各个问题的研究中，我通常都会采用不同的初值，发现不同初值会对应不同的迭代次数，另外针对问题一，我选用了不同的精度，发现迭代次数并没有很大的变化，因而一个初值的选定可以对迭代次数产生很大的影响，而精度的改变对迭代次数的影响很小。
对每个算法单独来看，显然选择初值不同对于迭代步数影响较大，对于找到根也会有影响。
因此应该先通过画图确定根的大致位置，给出在其附近的初值。
6、心得体会
在实现这三个算法的过程中，本身编程较易实现，最重要的是对算法本身的理解，只有真正理解算法的含义才能更快更好的实现程序。
、离散型最小二乘和连续型最小二乘问题
一、实验目的
掌握并能够利用离散型最小二乘和连续性最小二乘求解问题。
二、问题描述
1：
以函数
生成的6个数据点（0.25,23.1）,（1.0,1.68）,（1.5,1.0）,（2.0,0.84）,（2.4,0.826）,（5.0,1.2576）为例，运行程序得到不同阶对应的曲线拟合产生的多项式函数。
2.例题：
计算
f（x）=exp（x）在[-1，1]上的二、三次最佳平方逼近多项式。
并画图进行比较。
三、问题分析
问题1是离散最小二乘问题。
最小离散最小二乘就是根据一批有误的数据（
）i=1,2，…，n确定参数，并通过均方误差来比较曲线拟合的优劣，在本题中通过画图来比较不同阶方程拟合效果的优劣。
选择两种方法实现离散最小二乘。
方法一，建立normalequation（法方程组），求解k次多项式系数。
法方程组构造方法：
方法二：
由于在matlab中存在ployfit函数，可以即为方便的用k次多项式拟合。
问题2是一个连续型最小二乘法的应用实例。
对于给定的函数
，如果存在
使得
则称S*（x）是f（x）在集合
中的最佳平方逼近函数。
显然，求最佳平方逼近函数
的问题可归结为求它的系数
，使多元函数
取得极小值，也即点（
）是I（a0,…，an）的极点。
由于I（a0,a1,…，an）是关于a0,a1,…，an的二次函数，利用多元函数取得极值的必要条件，
（k=0,1,2,…,n）
即
得方程组
如采用函数内积记号
那么，方程组可以简写为
……………...
（1）
这是一个包含n+1个未知元a0,a1,…,an的n+1阶线性代数方程组，写成矩阵形式为
…………
（2）
此方程组叫做求aj（j=0,1,2,…,n）的法方程组。
显然，其系数行列式就是克莱姆行列式Gn=Gn（0,1,…,n）。
由于0,1,…,n线性无关，故Gn0，于是上述方程组存在唯一解
。
从而肯定了函数f（x）在
中如果存在最佳平方逼近函数，则必是
……………………………...（3）
四、实验程序
问题一
法一：
离散最小二乘法
functionf=normalequation（n）
symst
r=0;
xx=[0.25,1.0,1.5,2.0,2.4,5.0];
yy=[23.1,1.68,1.0,0.84,0.826,1.2576];x=zeros（n,n）;y=zeros（n,1）;
fori=1:
n
fork=1:
6
x（i,i）=xx（k）.^（2*i-2）+x（i,i）;
y（i,1）=yy（k）.*xx（k）.^（i-1）+y（i,1）;
end
end
fori=2:
n
forj=1:
i-1
fork=1:
6
x（i,j）=xx（k）.^（i+j-2）+x（i,j）;
x（j,i）=x（i,j）;
end
end
end
z=x\y;
fori=1:
n
r=z（i,1）*t^（i-1）+r;
end
vpa（r,5）
法二:
用matlab中的ployfit函数对k次多项式进行拟合。
建立number2.m文件，带入如下：
x=[0.25,1.0,1.5,2.0,2.4,5.0];
y=[23.1,1.68,1.0,0.84,0.826,1.2576];
plot（x,y,'o'）
p1=polyfit（x,y,1）%利用线性拟合
xi=-0.25:
0.01:
6.0;
y1=polyval（p1,xi）;
plot（x,y,'o',xi,y1,'r:
'）;
holdon;
p2=polyfit（x,y,2）%二次拟合
y2=polyval（p2,xi）;
plot（x,y,'o',xi,y2,'m'）;
holdon;
p3=polyfit（x,y,3）%三次拟合
y3=polyval（p3,xi）;
plot（x,y,'o',xi,y3,'b:
'）;
holdon;
p4=polyfit（x,y,4）%四次拟合
y4=polyval（p4,xi）;
plot（x,y,'o',xi,y4,'c'）;
holdon;
p5=polyfit（x,y,5）%五次拟合
y5=polyval（p5,xi）;
plot（x,y,'o',xi,y5,'g'）;
holdoff;
legend（'初始点','y=-2.3840*x+9.8499','初始点','y=2.1793*x^2-15.8845*x+22.3999','初始点','y=-2.0143*x^3+18.3499*x^2-45.2167*x+32.7386','初始点','y=1.4828*x^4-16.0435*x^3+56.3154*x^2-79.4994*x+39.6688','初始点','y=-0.9961*x^5+12.1743*x^4-55.2367*x^3+116.6262*x^2-116.7266*x+45.8090'）;
问题二：
最佳平方逼近算法
（1）输入被逼近函数f（x）和对应的逼近区间[a,b]并选择逼近函数系{
}和权函数；
（2）解方程组
（1）或
（2），其中方程组的系数矩阵和右端的项由式（3）得到；
（3）由式（3）得到函数的最佳平方逼近。
将上述算法编写成MATLAB程序共需三个程序：
第一个程序（函数名Sapproach.m）
计算最佳逼近函数的系数：
functionS=Sapproach（a,b,n）%定义逼近函数
globali;globalj;
ifnargin<3n=1;end%判断
X=zeros（n+1,n+1）;
fori=0:
n
forj=0:
n;
X（i+1,j+1）=quad（@fan,a,b）;%求fan积分
end
end
Y=zeros（n+1,1）;
fori=0:
n
Y（i+1）=quad（@yb,a,b）;%求yb积分
end
s=X\Y
第二个程序（函数名：
fan.m）
functiony=fan（x）
globali;globalj;
y=（poly（x,i））.*poly（x,j）;
;
第三个程序（函数名：
yb.m）
functiony=yb（x）
globali;
y=（poly（x,i））.*exp（x）;
编写多项式函数
functiony=poly（x,k）%多项式函数
ifk==0
y=ones（size（x））;
else
y=x.^k;
end
5、实验结果及图像
问题一：
（1）最小二乘法normalequation.m运行结果为：
线性拟合：
二次拟合：
三次拟合：
四次拟合：
五次拟合：
用方法二得到的结果如下：
number2.m运行结果为：
因此得到
线性拟合多项式为：
二次拟合多项式为：
三次拟合多项式为：
四次拟合多项式为：
五次拟合多项式为：
绘制的各多项式拟合图像为：
问题2：
当求二次逼近时，有以下结果：
绘制两者的图形：
>>fun='exp（x）';
>>fplot（fun,[-1,1]）
>>holdon
>>xi=-1:
0.1:
1;
>>yi=polyval（S,xi）;
>>plot（xi,yi,'r:
'）
得到以下结果：
当三次逼近时，有以下结果
绘制图形如下：
>>fun='exp（x）';
>>fplot（fun,[-1,1]）
>>holdon
>>xi=-1:
0.1:
1;
>>yi=polyval（S,xi）;
>>plot（xi,yi,'r:
'）
得到如下所得图形：
六、结果分析
显然离散最小二乘中次数较高拟合的效果较好，但由于本问题中初始点较少，并不能完全反应函数本身的特点，同时在本问题中，选择了两种方式，结果接近，也可以互相验证正确性。
在连续最小二乘法的实验中较顺利的达到了预期的结果。
从试验结果看出三次逼近没有二次逼近效果理想，验证了最佳平方逼近理论。
七、心得体会
在离散最小二乘与连续最小二乘程序设计的过程中，要么通过均差来表示函数拟合的优劣，要么通过图像来评价，本问题中选择了图像，更加清晰直观。
、数值积分
1、实验目的
熟悉并掌握数值积分的方法，重要训练复化梯形公式，复化simpson公式以及romberg积分。
2、问题描述
问题三数值积分椭圆周长的计算。
考虑椭圆
为计算其周长,只要计算其第一象限的长度即可.
用参数方程可以表示为
计算公式为
为计算方便,我们可以令
即计算下面的积分
（
可以归结为上面的形式）
采用复化梯形公式，复化Simpson公式以及Romberg积分的方法计算积分
给出通用程序，该通用程序可以计算任何一个函数在任意一个区间在给定的精度下的数值积分。
程序输出为计算出的数值积分值以及计算函数值的次数。
利用你的程序计算5个椭圆的周长。
3、算法介绍
首先利用给出的各迭代公式，设计程序。
在matlab对话框中输入要计算的函数，给出区间和精度。
复化梯形的迭代公式为：
复化simpson迭代公式为：
Romberg迭代公式为：
4、算法程序
复化梯形公式和复化simpson公式
我们放在jifen.m中。
（%标记后的程序可用来把b看为变量时的算法实现）
%复化梯形公式
functiony=jifen（f,n,a,b）（说明：
f表示任一函数，n精度，a，b为区间）
fi=f（a）+f（b）;
h=（b-a）/n;
d=1;
%functionf=jifen（n,a,b,c）
%symst
%y=sqrt（1+（c^2-1）*cos（t）^2）;
%ya=subs（y,t,a）;
%yb=subs（y,t,b）;
%fi=ya+yb;
fori=1:
n-1
x=a+i*h;
fi=fi+2*f（x）;
d=d+1;
%yx=subs（y,t,x）;
%fi=fi+2*yx;
end
f4=h/2*fi,d
%复化simposon公式
f1=0;
f2=0;
dd=1;
fori=1:
n-1
dd=dd+1;
ifrem（i,2）~=0;
x1=a+i*h;
f1=f1+f（x1）;
elserem（i,2）==0;
x2=a+i*h;
f2=f2+f（x2）;
end
end
f3=（h/3）*（f（a）+4*f1+2*f2+f（b））,dd
romberg积分
建立romberg.m文件
functiony=romberg（f,n,a,b）（说明：
f表示任一函数，n精度，a，b为区间）
z=zeros（n,n）;
h=b-a;
z（1,1）=（h/2）*（f（a）+f（b））;
f1=0;
fori=2:
n
fork=1:
2^（i-2）
f1=f1+f（a+（k-0.5）*h）;
end
z（i,1）=0.5*z（i-1,1）+0.5*h*f1;
h=h/2;
f1=0;
forj=2:
i
z（i,j）=z（i,j-1）+（z（i,j-1）-z（i-1,j-1））/（4^（j-1）-1）;
end
end
z,n
5、实验结果
（1）对于复化梯形公式和复化simpson公式，我们运行下列语句并得到结果：
题中将椭圆中的未知量a取为1，b取为0.5。
f4为复化梯形公式得到的1/4个椭圆周长，f3为复化simpson公式得到的1/4个椭圆周长。
从而得到椭圆的周长为4*1.2111=4.8444。
（2）对于romberg，运行下列语句并最终得到结果为：
题中将椭圆中的未知量
升级会员 