南京邮电大学数值计算实践报告Word下载.docx
- 文档编号:4592044
- 上传时间:2023-05-03
- 格式:DOCX
- 页数:47
- 大小:1.22MB
南京邮电大学数值计算实践报告Word下载.docx
《南京邮电大学数值计算实践报告Word下载.docx》由会员分享,可在线阅读,更多相关《南京邮电大学数值计算实践报告Word下载.docx(47页珍藏版)》请在冰点文库上搜索。
其中
是一阶均差和二阶均差。
收敛速度比割线法更接近于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
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);
建立n_df.m的源程序来方程的一阶导数的函数,源程序代码为:
functiony=n_df(a,n,x)%方程的一阶导数的函数
n
y=y+a(i)*(n+1-i)*x^(n-i);
在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
z=x-(feval(f,x)*(x-y))/(feval(f,x)-feval(f,y));
y=x;
x=z;
i
抛物线
建立paowuxian.m的源程序,源程序代码为:
functionx=paowuxian(f,x0,x1,x2,e)
x=x2;
y=x1;
z=x0;
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;
x=xi;
i
研究一:
只改变初值
由上述结果可知,方程的解在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软件中执行下列语句并得到的最终结果截图:
根据问题1中的抛物线法的程序,在matlab软件中执行下列语句并得到的最终结果截图:
、
研究结论
在只改变初值时,当初值定得越靠近初值,迭代次数就越少。
在只改变精度时,当精度越来越大时,迭代次数并几乎不变。
综上所述,初值对迭代次数的影响比较大,精度对迭代次数影响不大。
问题二
问题描述
问题分析
仍然利用
(1)中方法求解这一问题,并利用图解法找到初值,通过观察图像,将newton法初值设为:
0.1,割线法初值设为:
0,0.2。
抛物线法初值设为:
0,0.1,0.2。
图像见下图:
Newton法
问题一的牛顿法的求解只适用于线性方程,所以在问题二中用其他方法来求解方程。
建立newton1.m源程序,源程序代码为:
functionx=newton1(fn,dfn,x0,e)
ifnargin<
x=x0;
x0=x+2*e;
i=1;
whileabs(x0-x)>
i=i+1;
x0=x;
x=x0-feval(fn,x0)/feval(dfn,x0);
在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。
初值为0,0.3时
初值为0,0.1,0.2时
5、计算结果及分析
0.224
0.1466
0.224和
0.146577258557101
问题一
迭代步数
问题二
问题三
迭代步数
问题四
6
5
37
89
8
5
50
1310
9
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:
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);
fori=2:
forj=1:
i-1
x(i,j)=xx(k).^(i+j-2)+x(i,j);
x(j,i)=x(i,j);
z=x\y;
r=z(i,1)*t^(i-1)+r;
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);
xi,y1,'
r:
'
);
holdon;
p2=polyfit(x,y,2)%二次拟合
y2=polyval(p2,xi);
xi,y2,'
m'
p3=polyfit(x,y,3)%三次拟合
y3=polyval(p3,xi);
xi,y3,'
b:
p4=polyfit(x,y,4)%四次拟合
y4=polyval(p4,xi);
xi,y4,'
c'
p5=polyfit(x,y,5)%五次拟合
y5=polyval(p5,xi);
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;
3n=1;
end%判断
X=zeros(n+1,n+1);
fori=0:
forj=0:
n;
X(i+1,j+1)=quad(@fan,a,b);
%求fan积分
end
Y=zeros(n+1,1);
Y(i+1)=quad(@yb,a,b);
%求yb积分
end
s=X\Y
第二个程序(函数名:
fan.m)
functiony=fan(x)
y=(poly(x,i)).*poly(x,j);
;
第三个程序(函数名:
yb.m)
functiony=yb(x)
y=(poly(x,i)).*exp(x);
编写多项式函数
functiony=poly(x,k)%多项式函数
ifk==0
y=ones(size(x));
else
y=x.^k;
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,'
得到以下结果:
当三次逼近时,有以下结果
绘制图形如下:
得到如下所得图形:
六、结果分析
显然离散最小二乘中次数较高拟合的效果较好,但由于本问题中初始点较少,并不能完全反应函数本身的特点,同时在本问题中,选择了两种方式,结果接近,也可以互相验证正确性。
在连续最小二乘法的实验中较顺利的达到了预期的结果。
从试验结果看出三次逼近没有二次逼近效果理想,验证了最佳平方逼近理论。
七、心得体会
在离散最小二乘与连续最小二乘程序设计的过程中,要么通过均差来表示函数拟合的优劣,要么通过图像来评价,本问题中选择了图像,更加清晰直观。
、数值积分
熟悉并掌握数值积分的方法,重要训练复化梯形公式,复化simpson公式以及romberg积分。
问题三数值积分椭圆周长的计算。
考虑椭圆
为计算其周长,只要计算其第一象限的长度即可.
用参数方程可以表示为
计算公式为
为计算方便,我们可以令
即计算下面的积分
(
可以归结为上面的形式)
采用复化梯形公式,复化Simpson公式以及Romberg积分的方法计算积分
给出通用程序,该通用程序可以计算任何一个函数在任意一个区间在给定的精度下的数值积分。
程序输出为计算出的数值积分值以及计算函数值的次数。
利用你的程序计算5个椭圆的周长。
首先利用给出的各迭代公式,设计程序。
在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;
n-1
x=a+i*h;
fi=fi+2*f(x);
d=d+1;
%yx=subs(y,t,x);
%fi=fi+2*yx;
f4=h/2*fi,d
%复化simposon公式
f1=0;
f2=0;
dd=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);
f3=(h/3)*(f(a)+4*f1+2*f2+f(b)),dd
romberg积分
建立romberg.m文件
functiony=romberg(f,n,a,b)(说明:
z=zeros(n,n);
h=b-a;
z(1,1)=(h/2)*(f(a)+f(b));
2^(i-2)
f1=f1+f(a+(k-0.5)*h);
z(i,1)=0.5*z(i-1,1)+0.5*h*f1;
h=h/2;
f1=0;
forj=2:
z(i,j)=z(i,j-1)+(z(i,j-1)-z(i-1,j-1))/(4^(j-1)-1);
z,n
5、实验结果
(1)对于复化梯形公式和复化simpson公式,我们运行下列语句并得到结果:
题中将椭圆中的未知量a取为1,b取为0.5。
f4为复化梯形公式得到的1/4个椭圆周长,f3为复化simpson公式得到的1/4个椭圆周长。
从而得到椭圆的周长为4*1.2111=4.8444。
(2)对于romberg,运行下列语句并最终得到结果为:
题中将椭圆中的未知量a
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 南京 邮电大学 数值 计算 实践 报告
![提示](https://static.bingdoc.com/images/bang_tan.gif)