数值分析上机报告.docx
《数值分析上机报告.docx》由会员分享，可在线阅读，更多相关《数值分析上机报告.docx（23页珍藏版）》请在冰点文库上搜索。
数值分析上机报告
数值分析上机实习报告
专
业：
土木工程
班
级：
学
号：
姓
名：
指导老师：
联系电话：
2015.12.12
序言
随着本学期逐渐接近尾声，我也逐渐掌握了数值分析的一些基本理论•本次上机作业是理
论与实践的结合•
本次作业使用了matlab与C++两种语言•其中matlab具有编程效率高，用户使用方便，方便的绘图功能的优点。
而C++是一种基本的编程语言，在实际的工程中也有广泛的应用。
本次作业根据题目的特点，结合两种语言各自的优势，采用了不同的方法。
其中牛顿法，Steffensen加速法采用了c语言。
插值与多项式拟合使用了两种语言。
Runge-Kutta4阶算法
仅使用了matlab编程。
本次作业注重问题的计算过程，分析总结，及编程。
由于所涉及原理课本均有详细陈述，
在此不再赘述。
第一题3
1.1题目3
1.2计算过程和结果3
1.3结果分析3
第二题4
2.1题目4
2.2计算过程和结果4
2.3结果分析8
第三题8
3.1题目8
3.2问题求解及过程8
3.3结果分析9
总结10
附件11
第一题11
1.1.1第一问牛顿法11
1.1.2第一问牛顿-Steffensen法11
1.2.1第二问牛顿法12
1.2.2第二问牛顿-Steffensen法13
第二题14
2.1.1最小二乘法求解14
2.2.1拉格朗日差值多项式拟合15
2.2.2牛顿插值15
第三题17
3.1.1Runge-Kutta4阶算法17
第一题
1.1题目
分别用牛顿法，及基于牛顿算法下的Steffensen加速法
⑴求ln（x+sinx）=0的根。
初值x0分别取0.1,1,1.5,2,4进行计算。
（2）求sinx=0的根。
初值x0分别取1，1.4，1.6,1.8，3进行计算。
分析其中遇到的现象与问题。
1・2计算过程和结果
1.对方程In（x+sinx）=O，可求解x+sinx=1的解。
使用牛顿法，令f（x）=x•sinx-1，则
f（x）=「cosx，直至|xk..-xk|：
：
：
110-时，结束迭代；然后再使用基于牛顿法的Steffensen
加速法进行计算，直至|Xk1-Xk卜：
110占时，结束迭代。
其迭代结果与迭代次数如下表所
示（注N1为牛顿法迭代次数，N2为基于牛顿法Steffensen加速法迭代次数）：
x0
0.1
1
1.5
2
4
Newton迭代法
x
0.510973
0.510973
0.510973
0.510973
0.510973
N1
4
4
5
6
40
Newton—Steffensen
加速法
x
0.510973
0.510973
无收敛解
0.510973
0.510973
N2
3
3
4
8
2.对方程sinx=0，令f（x）=sinx，误差仍取|xk勺-xk卜：
110-5。
其迭代结果与迭代次数如
F表所示:
x0
1
1.4
1.6
1.8
3
Newton迭代法
x
0
3.14159
31.4159
6.28319
3.14159
N1
5
7
8
4
3
Newton—Steffensen力口速法
x
无收敛解
-3.14159
25.1327
6.28319
无收敛解
N2
4
6
3
1.3结果分析
1从牛顿与Steffensen加速法可以看出，牛顿一Steffensen加速法迭代速度明显快于牛顿加速法，从计算结果可以看出，x=0.510973处取得精确值，而在xO=O.1,1处收敛速度较快，说明牛顿法与基于牛顿法Steffensen加速法在单根附近有较快的收敛速度，迭代法是否收敛，与初始近似值xO的好坏很有关系。
2从sinx=0使用两种方法可以看出，Steffensen加速法迭代速度明显快于牛顿加速法，且
牛顿法和基于牛顿法Steffensen加速法在单根附近有较快的收敛速度，迭代法是否收敛，与
初始近似值xO的好坏很有关系，如在x0=1.6处，得到的收敛解很大，这是因为
f（1.6）=cos（1.6180/~）0，很难说明f（x）在此处是否发散或者是否收敛。
第二题
2.1题目
某过程测涉及两变量x和y,拟分别用插值多项式和多项式拟合给出其对应规律的近似多
项式，已知xi与yi之间的对应数据如下，xi=1,2,…,10
yi=34.658840.371914.6448-14.2721-13.357024.823475.2795
103.574397.484778.2392
（1）请用次数分别为3,4,5,6的多项式拟合并给出最好近似结果f（x）。
（2）请用插值多项式给出最好近似结果
下列数据为另外的对照记录，它们可以作为近似函数的评价参考数据。
xi=
Columns1through7
1.50001.9000
Columns8through14
2.3000
2.7000
3.1000
3.5000
3.9000
4.3000
4.7000
5.1000
5.5000
5.9000
6.3000
6.7000
Columns15through17
7.1000
7.5000
7.9000
Columns1through7
42.1498
41.4620
35.1182
24.3852
11.2732
-1.7813
-12.3006
Columns8through14
-18.1566
-17.9069-
■11.0226
2.0284
19.8549
40.3626
61.0840
Columns15through17
79.568893.7700102.3677
2.2计算过程和结果
第一问
编好拟合程序，计算结果如下：
（1）3次多项式拟合结果：
解得拟合多项式系数矩阵
e=[131.794399999964-94.47868142965119.333892249411-1.032633896659]T
故：
s3（x）=131.7944-94.4787X19.3338x2-1.0326x3
用MATLAB划出画出插值多项式函数曲线如下：
=3.2805e+003
「f（xj-S愉）]
2
i二
7.36797026998528
-0.3818X4:
用
（2）4次多项式拟合结果：
解得拟合多项式系数矩阵
e=[0.7449750004521073.53340190312746-42.14325642464905
-0.38184564393842]T
因此四次拟合多项式为s4（x）=0.7450-73.5334X-42.1433X2•7.3680x3
iw：
（3）5次多项式拟合结果：
解得拟合多项式系数矩阵
e=[-139.501853316416304.728173432415163.510704007400034.501972019191
3.0789000346780.098074705118]T
因此5次拟合多项式为：
S5（x）二-139.5019304.7282x-163.5107x234.5020x-3.0789x40.0981x5
（4）6次多项式拟合结果：
e=[-18.6990946534282966.37503804475497-0.86696651297073-16.89727867937389
5.11365302218162-0.540786038146500.01935941646399]T
因此6次拟合多项式为：
&（x）=-18.699166.375X-0.86702-16.89733+5.1137：
4-0.540&5+0.01946
2222比较不同次数的多项式拟合结果，可见2
说明6次多项式的最小二乘拟合的拟合程度最好•因此最好的拟合多项式为：
S6（x）--18.699166.3750x-0.8670x2-16.8973x35.1137x-0.5408x50.0194x6
第二问
用不同的插值方法进行插值多项式的求解。
（1）拉格朗日插值
求得9次拉格朗日插值多项式为：
N9（x）=100（-0.252042000119551.135********544x-0.8646529666577x20.50113073368872x3
-0.2240802436669636x4.0593********5-0.00879208298798x60.00072285644857x7
-0.00003088405259x80.00000053525132x9）
用上面的插值多项式计算出被插点数据x（参考点）对应的y:
Y=
42.384041.494735.074224.360111.2792-1.7683-12.2977-18.1626
-17.9118-11.02102.033319.856540.358461.079479.570993.7788
102.3713
计算所得到的值和参考值之间的方差为：
D1宀I]-yi）2=0.0589
（2）牛顿插值
使用9次牛顿插值多项式可得到同样的结果：
23
N9（x）=100（-0.252042000119551.13500691817544x_0.8646529666577x0.50113073368872x
_0.2240802436669636x4.0593********5-0.00879208298798x60.00072285644857x7-
89
0.00003088405259x0.00000053525132x）
N9（x）即为最佳插值多项式近似结果。
2.3结果分析
*2*2*2*2
1比较不同次数的多项式拟合结果，可见f-S62：
：
：
f-S52f-S42彳f-S32
说明次数越高最小二乘拟合的拟合程度最好。
2牛顿插值多项式跟拉格朗日插值的结果完全一致，这说明多项式是唯一的，因为插值节
点有10个，九次插值需要得到10个常系数，因此方程有唯一解，所以插值多项式是唯一的。
第二题
3.1题目
5用Runge-Kutta4阶算法对初值问题y/=-20*y，y（0）=1按不同步长求解，用于观察稳定区
间的作用，推荐两种步长h=0.1,0.2。
注：
此方程的精确解为：
y=e-20x
3.2问题求解及过程
（1）给定初值x为0,选取h=0.1进行计算，并对迭代过程和迭代图像进行分析。
计算结果如下：
请输入a的值：
0
请输入b的值：
1
请输入步长h:
0.1
X[1]=
0.000000
y[1]=
1.000000
yi[1]=
1.000000
error[1]=
0.000000
X[2]=
0.100000
y[2]=
0.333333
yi[2]=
0.135335
error[2]=
0.197998
x[3]=
0.200000
y[3]=
0.111111
yi[3]=
0.018316
error[3]=
0.092795
x[4]=
0.300000
y[4]=
0.037037
yi[4]=
0.002479
error[4]=
0.034558
x[5]=
0.400000
y[5]=
0.012346
yi[5]=
0.000335
error[5]=
0.012010
x[6]=
0.500000
y[6]=
0.004115
yi[6]=
0.000045
error[6]=
0.004070
x[7]=
0.600000
y[7]=
0.001372
yi[7]=
0.000006
error[7]=
0.001366
x[8]=
0.700000
y[8]=
0.000457
yi[8]=
0.000001
error[8]=
0.000456
x[9]=
0.800000
y[9]=
0.000152
yi[9]=
0.000000
error[9]=
0.000152
x[10]=
0.900000
y[10]=
0.000051
yi[10]=
0.000000
error[10]=
=0.000051
h=0.1
h=0.2
（2）给定初值x为0,,分别选取h=0.2进行计算，并对迭代过程和迭代图像进行分析。
计算结果如下
请输入a的值：
0
请输入b的值：
2
请输入步长h:
0.2
x[1]=
0.000000
y[1]=1.000000
yi[1]=
1.000000
error[1]=0.000000
x[2]=
0.200000
y[2]=5.000000
yi[2]=
0.018316
error[2]=4.981684
x[3]=
0.400000
y[3]=25.000000
yi[3]=
0.000335
error[3]=24.999665
x[4]=
0.600000
y[4]=125.000000
yi[4]=
0.000006
error[4]=124.999994
x[5]=
0.800000
y[5]=625.000000
yi[5]=
0.000000
error[5]=625.000000
x[6]=
1.000000
y[6]=3125.000000
yi[6]=
0.000000
error[6]=3125.000000
x[7]=
1.200000
y[7]=15625.000000
yi[7]=
0.000000
error[7]=15625.000000
x[8]=
1.400000
y[8]=78125.000000
yi[8]=
0.000000
error[8]=78125.000000
x[9]=
1.600000
y[9]=390625.000000
yi[9]=
0.000000
error[9]=390625.000000
x[10]=
1.800000
y[10]=1953125.000000
yi[10]=
0.000000
error[10]=1953125.000000
由图象和数据中可以看出，当步长为0.2时，迭代是不收敛的。
3.3结果分析
（1）用标准的四阶Runge-Kutta算法计算常微分方程时，不同的步长算法的稳定性有影响。
不够稳定的步长下面的计算，误差会越来越大，结果失真严重。
（2）一般情况下，步长越小，标准的四阶Runge-Kutta算法的稳定性越高，精度也越
总结
通过这次上机练习，让我对数值分析一些方法及其理论有了更深层次的理解，通过对
比，使我认识到各种方法的优缺点，并使我认识到学习数值分析中的不足，时的修补了自己知识上的漏洞，以便在复习与以后的工作中得以加强。
同时，通过本次上机实习，复习了本科阶段所学习的C++编程，提高了编程的熟练度，并且学习了一门新的语言-matlab，这些对我以后解决实际遇到的工程问题以及数学问题开辟了一个新的道路。
第一题
1.1.1第一问牛顿法
#include
#include
#defineEPS1.0e-5
intmain（）
{
doublex=0.0;
doubley=0.0;
intn=0;
cout«"请输入初值X0"<cin>>x;
do
{
y=x;
x=x-（x+sin（x）-1）/（1+cos（x））;
n++;
}
while（fabs（x-y）>=EPS）;
cout«"X"v}
1.1.2第一问牛顿-Steffensen法
#include
#include
#defineEPS1.0e-5intmain（）
doublex=0.0;
doubley=0.0;
doublez=0.0;
doublea=0.0;
intn=0;
cout«"请输入初值x0"<cin>>x;
do
{
a=x;
y=x-（x+sin（x）-1）/（1+cos（x））;
z=y-（y+sin（y）-1）/（1+cos（y））;
x=x-（y-x）*（y-x）/（z-2*y+x）;
n++;
}
while（fabs（x-a）>=EPS）;cout<<"x"<}
1.2.1第二问牛顿法
#include
#include
#defineEPS1.0e-5
intmain（）
{
doublex=0.0;
doubley=0.0;
intn=0;
cout«"请输入初值X0"«endl;
cin>>x;
do
{
y=x；
x=x-tan（x）;
n++;
}
while（fabs（x-y）>=EPS）;
cout<<"X"vvn}
1.2.2第二问牛顿-Steffensen法
#include
#include
#defineEPS1.0e-6
intmain（）
{
doublex=0.0;
doubley=0.0;
doublez=0.0;
doublea=0.0;
intn=0;
cout«"请输入初值x0"<cin>>x;
do
{
a=x;
y=x-tan（x）;
z=y-tan（y）;
x=x-（y-x）*（y-x）/（z-2*y+x）;
n++;
}
while（fabs（x-a）>=EPS）;
cout«"x"vcout«"迭代次数N="vv*}
第二题
A:
拟合多项式各项常
1范数
2.1.1最小二乘法求解
function[Anorm]=min2mul（x,y,n,x_count,y_count）%x,y:
原始数据点、
数项系数，n：
拟合次数x_count、y_count拟合参考点norm:
残差的fori=1:
n+1
forj=1:
n+1
a（i,j）=0;
fork=1:
length（x）
a（i,j）=a（i,j）+x（k）A（i-1）*x（k）q-1）;
end
end
end
fori=1:
n+1
b（i）=0;
fork=1:
length（y）
b（i）=b（i）+x（k）A（i-1）*y（k）;
end
end
A=a\b';
fori=1:
length（xi）
y_count_cal（i）=0;
forj=1:
n+1
y_count_cal（i）=y_count_cal（i）+A（j）*x_count（i）A（j-1）;
end
end
fori=1:
length（x）
y_cal（i）=0;
forj=1:
n+1
y_cal（i）=y_cal（i）+A（j）*x（i）A（j-i）;
end
end
Y=[yy_count];
Y_cal=[y_caly_count_cal];
norm=sum（（Y-Y_cal）.A2）
2.2.1拉格朗日差值多项式拟合
C++
#include
#include
#include
usingnamespacestd;
intmain（）{
intn,i,j;
cout<<"输入插值的次数"<cin»n;
floatX[n+1],Y[n+1],x,y,h[n],l;
cout<<"输入函数表"<for（i=0;i<=n;i++）{
cin>>X[i]>>Y[i];
}
y值（方便写出插值多项式）
for（i=0;i<=n;i++）{//求拉格朗日I的分子项乘以相应的
h[i]=1;
for（j=0;j<=n;j++）{
if（j!
=i）
h[i]*=（X[i]-X[j]）;
}
cout<}
while（cin>>x）{//用插值近似函数计算
y=0;
for（i=0;i<=n;i++）{
l=1;
for（j=0;j<=n;j++）{
if（j!
=i）
l=l*（x-X[j]）/（X[i]-X[j]）;
}
y=y+l*Y[i];
}
cout<<"y="<}
return0;
}
2.2.2牛顿插值
#include
#include
usingnamespacestd;
#include
intmain（）{
intn,i,j;
floath,x,y,xO,t,temp1,temp2,T;
cout<<"输入插值节点个数"<floatX[n],Y[n];
cout<<"输入函数表"<cin>>X[i]>>Y[i];
}
h=X[1]-X[0];//步长while（cin>>x）{
cout<<"输入插值节点xO:
"<>x0;
t=（x-xO）/h;
y=0;for（i=0;ifor（j=i+1;jtemp2=
升级会员 