东南大学数值分析上机作业汇总.docx
《东南大学数值分析上机作业汇总.docx》由会员分享，可在线阅读，更多相关《东南大学数值分析上机作业汇总.docx（17页珍藏版）》请在冰点文库上搜索。
东南大学数值分析上机作业汇总
数值分析
上机报告
院系：
学号：
姓名：
作业1、舍入误差与有效数
设
，其精确值为
.
（1）编制按从小到大的顺序
，计算
的通用程序；
（2）编制按从大到小的顺序
，计算
的通用程序；
（3）按两种顺序分别计算
并指出有效位数；
（4）通过本上机你明白了什么？
程序：
1、函数文件cxdd.m
functionS=cxdd（N）
S=0;
i=2.0;
while（i<=N）
S=S+1.0/（i*i-1）;
i=i+1;
end
script运行结果（省略>>）：
S=cxdd（80）
S=
0.737577
2、函数文件cddx.m
functionS=cddx（N）
S=0;
fori=N:
-1:
2
S=S+1/（i*i-1）;
end
script运行结果（省略>>）：
S=cddx（80）
S=
0.737577
3、两种方法有效位数对比
精确值函数：
functionS=jqz（N）
S=0.5*（1.5-1.0/N-1.0/（N+1））;
script运行结果（省略>>）
N
S精确值
从小到大
从大到小
值
有效位数
值
有效位数
100
0.740050
0.740050
6
0.740049
6
10000
0.749900
0.749900
4
0.749852
4
1000000
0.749999
0.749999
6
0.749852
3
4、心得
本题重点体现了数值计算中“大数吃小数”的问题，由于计算机计算的截断特点，从大到小的计算会导致小数的有效数被忽略掉。
从题中可以看出，看出按不同的顺序计算的结果是不相同的，按从小到大的顺序计算的值与精确值吻合，而按从大到小的顺序计算的值与精确值有较大的误差。
计算机在进行数值计算时会出现“大数吃小数”的现象，导致计算结果的精度有所降低。
作业2、Newton迭代法
（1）给定初值x0及容许误差ε，编制Newton法解方程f（x）=0根的通用程序。
（2）给定方程f（x）=x3/3-x=0,易知其有三个根x1※=
，x2※=0，x3※=
。
由Newton方法的局部收敛性可知存在
＞0，当x0∈（
，
），Newton迭代序列收敛于根x2※，试确定尽可能大的
；
②试取若干个初始值，观察当x0∈（-∞，-1），（-1，
），（
，
），（
，1），（1，+∞）时，Newton序列是否收敛以及收敛于哪一个根。
（3）通过本上机题，你明白了什么？
1、通用程序函数文件
定义f（x）函数
functionf=fun（x）
f=x^3/3-x;
end
定义f（x）导函数
functionf=dfun（x）
f=x*x-1;
end
定义求近似解函数
function[f,n]=newton（x0,ep）
flag=1;
n=0;
while（flag==1）
x1=x0-fun（x0）/dfun（x0）;
n=n+1;
if（abs（x1-x0）<=ep||n>100000）
flag=0;
end
x0=x1;
end
f=x1;
end
script运行结果
clear;
x0=input（'请输入初始值x0：
'）;
ep=input（'请输入容许误差：
'）;
[f,n]=newton（x0,ep）;
fprintf（'方程的一个近似解为:
%f\n',x1）;
2、局部收敛性
（1）最大δ值文件
flag=1;
k=1;
x0=0;
whileflag==1
sigma=k*10^-6;
x0=sigma;
k=k+1;
m=0;
flag1=1;
whileflag1==1&&m<=10^3
x1=x0-fun（x0）/dfun（x0）;
ifabs（x1-x0）<10^-6
flag1=0;
end
m=m+1;
x0=x1;
end
if（flag1==1||abs（x0）>=10^-6）
flag=0;
end
end
fprintf（'最大值为：
%f\n',sigma）;
运行结果为：
最大值为：
0.774597
即得最大的δ为0.774597，Newton迭代序列收敛于根
=0的最大区间为（-0.774597，0.774597）。
（2）验证局部收敛性
●在x0∈（-∞，-1）区间，取以下初值，分别调用newton.m函数文件，得到结果如下：
X0
X1
迭代次数
-100
-1.732051
15
-20
-1.732051
11
-5
-1.732051
8
-1.5
-1.732051
5
结果显示，以上初值迭代序列均收敛于-1.732051，即根
。
显然，迭代格式初值的选择对于迭代的收敛速度是至关重要的，当初值接近真实值的时候，迭代次数减少。
●在x0∈（-1，
）区间，取以下初值，分别调用newton.m函数文件，得到结果如下：
X0
X1
迭代次数
-0.95
1.732051
9
-0.85
1.732051
6
-0.80
1.732051
10
-0.78
1.732051
15
计算结果显示，迭代序列局部收敛于1.730251，即根
。
●在x0∈（
，
）区间，取以下初值，分别调用newton.m函数文件，得到结果如下：
X0
X1
迭代次数
-0.70
0.000000
5
-0.20
0.000000
3
-0.05
0.000000
3
0.05
0.000000
3
0.20
0.000000
3
0.70
0.000000
5
由newton1.m的运行过程表明，在整个区间上均收敛于0，即根
。
●在x0∈（
，1）区间，取以下初值，分别调用newton.m函数文件，得到结果如下：
X0
X1
迭代次数
0.80
-1.732051
10
0.90
-1.732051
7
0.95
-1.732051
9
0.98
-1.732051
12
计算结果显示，迭代序列局部收敛于-1.732051，即根
。
●在x0∈（1，+∞）区间，取以下初值，分别调用newton.m函数文件，得到结果如下：
X0
X1
迭代次数
1.5
1.732051
5
5
1.732051
8
20
1.732051
11
100
1.732051
15
结果显示，以上初值迭代序列均收敛于1.732051，即根
。
综上所述：
（-∞，-1）区间收敛于-1.73205，（-1，δ）区间局部收敛于1.73205，局部收敛于-1.73205，（-δ，δ）区间收敛于0，（δ，1）区间类似于（-1，δ）区间，（1,∞）收敛于1.73205。
3、心得
牛顿迭代法对于初值的选择要求较高，因此，在牛顿迭代时可现通过简单迭代法寻找相对准确一些的值来进行牛顿迭代。
对于方程有多解的问题，Newton法求方程根时，牛顿迭代要考虑局部收敛的问题，迭代序列收敛于某一个根有一定的区间限制，在一个区间上，可能会局部收敛于不同的根。
作业3、列主元素Gauss消去法
对于某电路的分析，归结为求解线性方程组RI=V。
32-13000-10000
-1335-90-110000
0-931-1000000
R=000-3057-70-50
0000-747-3000
00000-304100
0000-50027-2
000-9000-229
VT=[-15，27，-23，0，-20,12，-7,7,10]T
（1）编制解n阶线性方程组Ax=b的列主元Gauss消去法的通用程序；
（2）用所编程序解线性方程组RI=V，并打印出解向量，保留5位有效数字；
（3）在本编程之中，你提高了那些编程能力。
1、列主元Gauss消去法的通用程序
函数：
找每列的主元的函数
functionB=zhuyuan（B,t,N,M）
fori=0:
N-1-t
ifB（N-i,t）>B（N-i-1,t）
c=zeros（1,M）;
forj=1:
M
c（j）=B（N-i,j）;
B（N-i,j）=B（N-i-1,j）;
B（N-i-1,j）=c（j）;
end
end
end
进行列消去的函数
functionB=xiaoqu（B,t,N,M）
fori=t+1:
N
l=B（i,t）/B（t,t）;
forj=t:
M
B（i,j）=B（i,j）-l*B（t,j）;
end
end
进行三角矩阵下的解函数
functionX=jie（X,B,N,M）
fori=1:
N-1
s=B（N-i,M）;
forj=N-i+1:
N
s=s-B（N-i,j）*X（j）;
end
X（N-i）=s/B（N-i,N-i）;
end
执行主程序：
N=input（'请输入线性方程组的阶数:
N='）;
M=input（'请输入增广矩阵阶数:
M='）;
b=zeros（1,N）;
A=zeros（N,N）;
A=input（'请输入系数矩阵:
'）;
b（1,:
）=input（'请输入方程组的右端向量:
'）;
b=b';
B=[A,b];
fort=1:
N-1B=zhuyuan（B,t,N,M）;
B=xiaoqu（B,t,N,M）;
end
X=zeros（N,1）;
X（N）=B（N,M）/B（N,N）;
X=jie（X,B,N,M）;
X
2、解题中线性方程组
执行程序，输入矩阵A（即题中的矩阵R）和列向量b（即题中的V）,如下：
请输入线性方程组的阶数:
n=9
请输入增广矩阵阶数:
m=10
请输入系数矩阵A：
A=[31,-13,0,0,0,-10,0,0,0,-15;-13,35,-9,0,-11,0,0,0,0,27;0,-9,31,-10,0,0,0,0,0,-23;0,0,-10,79,-30,0,0,0,-9,0;0,0,0,-30,57,-7,0,-5,0,-20;0,0,0,0,-7,47,-30,0,0,12;0,0,0,0,0,-30,41,0,0,-7;0,0,0,0,-5,0,0,27,-2,7;0,0,0,-9,0,0,0,-2,29,10];
请输入方程组的右端向量b：
[-1527-230-2012-7710]
得到如下结果：
A=
31-13000-10000-15
-1335-90-11000027
0-931-1000000-23
00-1079-30000-90
000-3057-70-50-20
0000-747-300012
00000-304100-7
0000-50027-27
000-9000-22910
x=
-0.28923
0.34544
-0.71281
-0.22061
-0.43040
0.15431
-0.057823
0.20105
0.29023
3、心得
列主元Gauss最重要的就是如何通过找到最大主元，并交换行，如何进行消去，这需要很细心的循环，和很复杂的嵌套，通过本题的编程，感觉个人对于这种多层嵌套循环的处理能力提高了很多。
通过对软件的编程，也更加理解了Gauss消去法的实质。
作业4、三次样条插值函数
（1）编制求第一型3次样条插值函数的通用程序；
（2）已知汽车门曲线型值点的数据如下：
i
0
1
2
3
4
5
6
7
8
9
10
xi
0
1
2
3
4
5
6
7
8
9
10
yi
2.51
3.3
4.04
4.7
5.22
5.54
5.78
5.4
5.57
5.7
5.8
端点条件为
=0.8，
=0.2，用所编程序求车门的3次样条差值函数S（x），并打印出S（i+0.5），i=0,1······9。
1、第一型三次样条插值函数通用程序：
n=input（'请输入节点数n:
'）;
n=n+1;
xn=zeros（1,n）;
yn=zeros（1,n）;
xn（1,:
）=input（'请输入X的值：
'）;
yn（1,:
）=input（'请输入Y的值：
'）;
dy0=input（'请输入边界条件y（0）:
'）;
dyn=input（'请输入边界条件y（n）:
'）;
d=zeros（n,1）;
h=zeros（1,n-1）;
f1=zeros（1,n-1）;
f2=zeros（1,n-2）;
fori=1:
n-1
h（i）=xn（i+1）-xn（i）;%求一阶差商
f1（i）=（yn（i+1）-yn（i））/h（i）;
end
fori=2:
n-1
f2（i）=（f1（i）-f1（i-1））/（xn（i+1）-xn（i-1））;%求二阶差商
d（i）=6*f2（i）;
end
d
（1）=6*（f1
（1）-dy0）/h
（1）;
d（n）=6*（dyn-f1（n-1））/h（n-1）;
A=zeros（n）;%求M值
u=zeros（1,n-2）;
r=zeros（1,n-2）;
fori=1:
n-2
u（i）=h（i）/（h（i）+h（i+1））;
r（i）=1-u（i）;
end
A（1,2）=1;
A（n,n-1）=1;
fori=1:
n
A（i,i）=2;
end
fori=2:
n-1
A（i,i-1）=u（i-1）;
A（i,i+1）=r（i-1）;
end
M=A\d;
symsx
fori=1:
n-1%求节点插值
Sx（i）=collect（yn（i）+（f1（i）-（M（i）/3+M（i+1）/6）*h（i））*（x-xn（i））
+M（i）/2*（x-xn（i））^2+（M（i+1）-M（i））/（6*h（i））*（x-xn（i））^3）;
Sx（i）=vpa（Sx（i）,4）;
end
S=zeros（1,n-1）;
fori=1:
n-1
x=xn（i）+0.5;
S（i）=yn（i）+（f1（i）-（M（i）/3+M（i+1）/6）*h（i））*（x-xn（i））
+M（i）/2*（x-xn（i））^2+（M（i+1）-M（i））/（6*h（i））*（x-xn（i））^3;
end
disp（'S（x）='）;%结果输出
fori=1:
n-1
fprintf（'%s（%d,%d）\n',char（Sx（i））,xn（i）,xn（i+1））;
disp（'··············································'）;
end
disp（'S（i+0.5）'）;
disp（'ix（i+0.5）S（i+0.5）'）
fori=1:
n-1
fprintf（'%d%.4f%.4f\n',i,xn（i）+0.5,S（i））
end
2、数据输入及计算结果
请输入节点数n:
10
请输入X的值：
[012345678910]
请输入Y的值：
[2.513.304.044.705.225.545.785.405.575.705.80]
请输入边界条件y（0）:
0.8
请输入边界条件y（n）:
0.2
S（x）=
-0.008514*x^3-0.001486*x^2+0.8*x+2.51（0,1）
·············································
-0.0044579*x^3-0.013654*x^2+0.81217*x+2.5059（1,2）
·············································
-0.0036544*x^3-0.018475*x^2+0.82181*x+2.4995（2,3）
·············································
-0.040924*x^3+0.31696*x^2-0.18448*x+3.5058（3,4）
·············································
0.10735*x^3-1.4624*x^2+6.9328*x-5.9839（4,5）
·············································
-0.26848*x^3+4.1752*x^2-21.255*x+40.996（5,6）
·············································
0.42659*x^3-8.3361*x^2+53.813*x-109.14（6,7）
·············································
-0.26786*x^3+6.2474*x^2-48.272*x+129.06（7,8）
·············································
0.054872*x^3-1.4983*x^2+13.694*x-36.184（8,9）
·············································
0.058376*x^3-1.5929*x^2+14.545*x-38.738（9,10）
·············································
S（i+0.5）
ix（i+0.5）S（i+0.5）
10.50002.9086
21.50003.6784
32.50004.3815
43.50004.9882
54.50005.3833
65.50005.7237
76.50005.5944
87.50005.4299
98.50005.6598
109.50005.7323
升级会员 