数值计算课程设计.docx
《数值计算课程设计.docx》由会员分享，可在线阅读，更多相关《数值计算课程设计.docx（20页珍藏版）》请在冰点文库上搜索。
数值计算课程设计
数值计算课程设计报告
系别
信息与计算科学
专业
学号
姓名
指导教师
成绩
教师评语：
指导教师签字：
2011年7月13日
1绪论
数值分析是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。
MATLAB是主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。
它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案。
2用MATLAB进行插值计算
题目
下列数据点的插值
X
0
1
4
9
16
25
36
49
64
Y
0
1
2
3
4
5
6
7
8
可以得到平方根函数的近似，在区间[0,64]上作图。
（1）用这9个点作8次多项式插值
.
（2）用三次样条（第一边界条件）程序求
.
从得到的结果看在[0,64]上，哪个插值更精确；在区间[0,1]上，两种插值哪个更精确？
2.1程序文件
建立新的m文件unit2.m：
代码如下：
clc;clear;
x=[0,1,4,9,16,25,36,49,64];
y=[0:
8];
xi=[0:
64];
p=polyfit（x,y,8）;
yi=polyval（p,xi）;
subplot（2,2,1）;
plot（x,y,'x',xi,yi,'k'）
xlabel（'x轴'）;ylabel（'y轴'）;
legend（'插值点','8次差值多项式'）
yi=interp1（x,y,xi,'cubic'）;
subplot（2,2,2）;
plot（x,y,'x',xi,yi,'k'）
xlabel（'x轴'）;ylabel（'y轴'）;
legend（'插值点','三次样条插值'）
xi=[0:
0.01:
1];
p=polyfit（x,y,8）;
yi=polyval（p,xi）;
subplot（2,2,3）;
plot（xi,yi,'k'）
xlabel（'x轴'）;ylabel（'y轴'）;
legend（'8次差值多项式'）
yi=interp1（x,y,xi,'cubic'）;
subplot（2,2,4）;
plot（xi,yi,'k'）
xlabel（'x轴'）;ylabel（'y轴'）;
legend（'三次样条插值'）
图1
2.2图样
运行文件得到图1
2.3分析
比较两个函数的图像可知，在区间[0,64]上三次样条插值函数要更精确一些，在区间[0,1]上拉格朗日插值函数仍然不如三次样条插值函数更精确。
3多项式曲线拟合
题目
由实验给出数据表
X
0.0
0.1
0.2
0.3
0.5
0.8
1.0
Y
1.0
0.41
0.50
0.61
0.91
2.02
2.46
试求3次、4次多项式的曲线拟合，再根据曲线形状，求一个另外函数的拟合曲线，用图示数据曲线及相应的三种拟合曲线。
3.1程序文件
建立新的m文件unit3.m：
代码如下：
clc;clear;
x=[0.0,0.1,0.2,0.3,0.5,0.8,1.0];
y=[1.0,0.41,0.50,0.61,0.91,2.02,2.46];
a=polyfit（x,y,3）;
b=polyfit（x,y,4）;
xi=[0.0:
0.01:
1.0];
aa=polyval（a,xi）;
bb=polyval（b,xi）;
subplot（1,2,1）;
plot（x,y,'x',xi,aa,'k'）
xlabel（'x轴'）;ylabel（'y轴'）;
legend（'插值点','3次曲线拟合'）
subplot（1,2,2）;
plot（x,y,'x',xi,bb,'k'）
xlabel（'x轴'）;ylabel（'y轴'）;
legend（'插值点','4次曲线拟合'）
poly2str（a,'x'）
poly2str（b,'x'）
3.2图形和结果
运行程序结果如下：
ans=-6.6221x^3+12.8147x^2-4.6591x+0.92659
ans=2.8853x^4-12.3348x^3+16.2747x^2-5.2987x+0.94272
图像如下：
图2
3.25次多项式曲线拟合
下面进行的：
代码如下：
clc;clear;
x=[0.0,0.1,0.2,0.3,0.5,0.8,1.0];
y=[1.0,0.41,0.50,0.61,0.91,2.02,2.46];
a=polyfit（x,y,5）;
xi=[0.0:
0.01:
1.0];
aa=polyval（a,xi）;
plot（x,y,'x',xi,aa,'k'）
xlabel（'x轴'）;ylabel（'y轴'）;
legend（'插值点','5次曲线拟合'）
poly2str（a,'x'）运行程序结果如下：
ans=-79.3261x^5+195.4554x^4-172.7104x^3+69.0498x^2-11.0044x+0.99547
图像如下：
图3
由图像可知，次数越高，拟合越准确。
4数值积分
题目
用不同数值方法计算积分
.
（1）取不同的步长h.分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能呢个再被改善？
（2）用龙贝格求积计算完成问题
（1）.
（3）用自适应辛普森积分，使其精度达到
4.1复合梯形公式代码
建立新的m文件unit41.m：
代码如下：
clc,clear;
x=0.0000001:
0.1:
1;
y=x.^（1/2）.*log（x）;
Fx=trapz（x,y）
error=Fx+4/9
运行文件结果为：
Fx=-0.4123
error=0.0321
取不同步长h:
h=0.01，Fx=-0.4431，error=-0.0014；
h=0.001，Fx=-0.4444，error=-5.4875e-005；
h=0.0001，Fx=-0.4444，error=-2.0338e-006；
h=0.00001，Fx=-0.4444，error=-6.4496e-008.
可见，没有一个最小的h使精度不能再改变。
4.2复合辛普森公式代码
建立新的m文件unit42.m：
代码如下：
Fx=quad（inline（'x.^（1/2）.*log（x）'）,0.000001,1，0.0001）
error=Fx+4/9
运行文件结果为：
Fx=-0.4440
error=4.3273e-004
取不同的步长h：
h=0.00001，Fx=-0.4444，error=-6.3537e-005；
h=0.000001，Fx=-0.4444，error=-2.9938e-006；
h=0.0000001，Fx=-0.4444，error=-3.0657e-007；
h=0.0000000000001，Fx=-0.4444，error=-9.6548e-009；
h=0.00000000000001，Fx=-0.4444，error=-9.6548e-009。
可见，存在一个最小的h使精度不能在改善，且h在0.0000000000001与0.00000000000001之间。
4.3龙贝格算法
建立新的m文件unit43.m：
代码如下：
function[q,n]=unit43（f,a,b）;
M=1;
abs0=10;
k=0;
T=zeros（1,1）;
h=b-a;
T（1,1）=（h/2）*（subs（f,a）+subs（f,b））;
whileabs0>0.0001
k=k+1;
h=h/2;
p=0;
fori=1:
M
x=a+h*（2*i-1）;
p=p+subs（f,x）;
end
T（k+1,1）=T（k,1）/2+h*p;
M=2*M;
forj=1:
k
T（k+1,j+1）=（（4^j）*T（k+1,j）-T（k,j））/（4^j-1）;
end
abs0=abs（T（k+1,j+1）-T（k,j））;
end
q=T（k+1,k+1）;
n=k;
在命令窗口输入：
[Fx,n]=unit43（'sqrt（x）*log（x）',10^（-8）,1）
输出结果为：
Fx=-0.4444
n=9
4.4自适应辛普森算法
Fx=quad（inline（'x.^（1/2）.*log（x）'）,0.000001,1）
输出结果为：
Fx=-0.4444。
5解线性方程
题目
线性方程组Ax=b的A及b为
A=
，b=
，
则解x=
.用MATLAB内部函数求detA及A的所有特征值和cond（A）
.若令
A+δA=
，
求解（A+δA）（x+δx）=b,输出向量δx和
.从理论结果和实际计算两方面分析线性方程组Ax=b解的相对误差
及A的相对误差
的关系.
5.1求detA及A的所有特征值和cond（A）
输入程序：
A=[10787;7565;86109;75910]
det（A）
[V,D]=eig（A）
norm（A,2）
输出结果为：
A=
10787
7565
86109
75910
ans=1
V=
0.5016-0.30170.61490.5286
-0.83040.09330.39630.3803
0.20860.7603-0.27160.5520
-0.1237-0.5676-0.62540.5209
D=
0.0102000
00.843100
003.85810
00030.2887
ans=30.2887
则得到detA=1
A的特征值为0.01010.84313.858130.2887
=30.2887
6解线性方程组的迭代法
题目
给出线性方程组
其中系数矩阵
为希尔伯矩阵：
，
，i=1,2……n.
假设
，
，若取n=6,8,10，分别用雅克比迭代及SOR迭代（w=1,1.25,1.5）求解.比较计算结果.
6.1用雅可比迭代法求解
建立新的m文件unit61.m：
代码如下：
n=6;
H=hilb（n）;
b=H*ones（n,1）;
e=0.00001;
fori=1:
n
ifH（i,i）==0'对角元为零，不能求解'
return
end
end
x=zeros（n,1）;
k=0;
kend=10000;
r=1;
whilek<=kend&r>e
x0=x;
fori=1:
n
s=0;
forj=1:
i-1s=s+H（i,j）*x0（j）;
end
forj=i+1:
ns=s+H（i,j）*x0（j）;
end
x（i）=b（i）/H（i,i）-s/H（i,i）;
end
r=norm（x-x0,inf）;
k=k+1;
end
ifk>kend'迭代不收敛，失败'
else'求解成功'
x
k
end
输出结果为：
ans=迭代不收敛，失败
6.2用SOR迭代法求解
建立新的m文件unit62.m：
代码如下：
functions=unit62（n,w）;
H=hilb（n）;
b=H*ones（n,1）;
e=0.00001;
fori=1:
n
ifH（i,i）==0
'对角元为零，不能求解'
return
end
end
x=zeros（n,1）;
k=0;
kend=10000;
r=1;
whilek<=kend&r>e
x0=x;
fori=1:
n
s=0;
forj=1:
i-1s=s+H（i,j）*x（j）;
end
forj=i+1:
ns=s+H（i,j）*x0（j）;
end
x（i）=（1-w）*x0（i）+w/H（i,i）*（b（i）-s）;
end
r=norm（x-x0,inf）;
k=k+1;
end
ifk>kend'迭代不收敛，失败'
else'求解成功'
x
k
end
输入：
unit62（6,1）
输出结果为：
ans=求解成功
x=0.99931.01310.95371.03741.02960.9662
k=997
ans=0.6487
当n=6，w=1.25时，ans=求解成功，x=[0.9992,1.0131,0.9558,1.0375,1.0197,0.9741]，k=3121，ans=0.6480；
当n=6，w=1.5时，ans=求解成功，x=[0.9995,1.0081,0.9726,1.0232,1.0123,0.9839]，k=5635，ans=0.6471；
当n=8，w=1时，ans=求解成功，x=[0.9996,1.0039,0.9938,0.9824,1.0071，1.0200，1.0093,0.9790]，k=3808，ans=0.6601；
当n=8，w=1.25时，ans=求解成功x=[0.9998,1.0004,1.0103,0.9685,1.0109,1.0201,1.0100,0.9795]，k=3949，ans=0.6601；
当n=8，w=1.5时，ans=求解成功，x=[1.0000,0.9975,1.0231,0.9469,1.0275,1.0118,1.0158,0.9771]，k=4031，ans=0.6602；
当n=10，w=1时，ans=求解成功，x=[1.0001,0.9952,1.0275,0.9677,0.9816,1.0093,1.0241,1.0208,1.0018,0.9713]，k=2673，ans=0.6677,；
当n=10，w=1.25时，ans=求解成功，x=[1.0004,0.9904,1.0466,0.9421,0.9886,1.0123,1.0272,1.0214,1.0010,0.9694]，k=2562，ans=0.6677；
当n=10，w=1.5时，ans=求解成功，x=[1.0006,0.9857,1.0688,0.9028,1.0177,0.9993,1.0377,1.0169,1.0030,0.9669]，k=2463，ans=0.6679
7用MATLAB求解非线性方程的根
题目
求下列方程的实根
（1）
;
（2）
.
要求：
（1）设计一种不动点迭代法，要使迭代序列收敛，然后再用斯特芬森加速迭代，计算到
为止.
（2）用牛顿迭代，同样计算到
.输出迭代初值及各次迭代值和迭代次数k，比较方法的优劣.
7.1不动点迭代法
建立新的m文件unit71.m：
代码如下：
function[y,n]=unit71（x,eps）
ifnargin==1
eps=1.0e-8;
elseifnargin<1
error
return
end
x1=gg（x）;n=1;
while（norm（x1-x）>=1e-8）&&（n<=10000）
x=x1;x1=gg（x）;n=n+1;
end
y=x;
若是functionf=gg（x）f
（1）=（exp（x）-2-x^2）/（-3）;
在在命令窗口输入：
unit71
（1）
输出结果为：
n=15
ans=0.2575
若是functionf=gg（x）f
（1）=-（x^3+2*x^2-20）/10;
在命令窗口输入：
unit71
（1）
输出结果为：
n=10001
ans=0.5489
7.2斯特芬森加速迭代法
建立新的m文件unit72.m：
代码如下：
functiony=unit72（x0,eps）
ifnargin<2
eps=1.0e-8;
end
a=g（x0）;b=g（a）;
x1=x0-（a-x0）^2/（b-2*a+x0）;n=1;
while（abs（x1-x0）>=eps）&（n<=10000）
x0=x1;a=g（x0）;b=g（a）;
x1=x0-（a-x0）^2/（b-2*a+x0）;n=n+1;
end
x1
n
若是functiony=g（x）;y=（exp（x）-2-x^2）/（-3）;
在命令窗口输入：
unit72
（1）
输出结果为：
x1=0.2575
n=4
若是functiony=g（x）;y=-（x^3+2*x^2-20）/10
在命令窗口输入：
unit72
（1）
输出结果为：
x1=1.3688
n=5
7.3牛顿迭代法
建立新的m文件unit73.m：
代码如下：
functiony=unit73（x0,eps）
ifnargin<2
eps=1.0e-8;
end
x1=x0-fa（x0）/faa（x0）;n=1;
while（abs（x1-x0）>=eps）&（n<=10000）
x0=x1;x1=x0-fa（x0）/faa（x0）;n=n+1;
end
x1
n
若是functiony=fa（x）;y=x^2-3+2*x-exp（x）;
functiony=faa（x）;y=2*x-3-exp（x）;
在命令窗口输入：
unit73（0）
输出结果为：
x1=-3.0123
n=33
若是functiony=fa（x）;y=x^3+2*x^2+10*x-20;
functiony=faa（x）;y=3*x^2+4*x+10;
在命令窗口输入：
unit73（0）
输出结果为：
x1=1.3688
n=6
通过三种方法的比较，可以看出斯特芬森加速迭代法迭代次数最少，方法最优，牛顿迭代法次之，而不动点迭代法次数最多，不动点迭代法针对收敛函数求解，所以在不动点迭代法求解之前要先判断函数敛散性。
牛顿迭代法因其效率高而应用最为广泛，但是它对重很收敛很慢，而且对初值要求严格。
结论
数值分析是研究如何用计算机解决实际问题的课程．将MATLAB与数值分析课程结合起来，开阔了思路，拓展了解决问题的方法，取得了较好的效果．MATLAB在数值计算中发挥着十分重要的作用．
通过这段时间对MATLAB的学习，我对MATLAB有了更进一步的了解，能够运用MATLAB来编写简单的函数和程序，也能进行符号计算，尤其是我可以通过MATLAB来解决数值分析的一些问题．我对计算机解决实际问题有了一个全新的认识．我会继续努力地学习MATLAB软件，掌握更多的相关知识．
参考文献
[1]李庆杨,王能超,易大义.数值分析（第5版）[M].北京:
清华大学出版社,2006.
[2]胡良剑,孙晓君,MATLAB数学实验[M].北京:
高等教育出版社,2006.
升级会员 