数值分析MATLAB上机实验.docx
《数值分析MATLAB上机实验.docx》由会员分享，可在线阅读，更多相关《数值分析MATLAB上机实验.docx（29页珍藏版）》请在冰点文库上搜索。
数值分析MATLAB上机实验
数值分析实习报告
姓名：
gestepoA
学号：
201*******
班级：
***班
序言
随着计算机技术的迅速发展，数值分析在工程技术领域中的应用越来越广泛，并且成为数学与计算机之间的桥梁。
要解决工程问题，往往需要处理很多数学模型，不仅要研究各种数学问题的数值解法，同时也要分析所用的数值解法在理论上的合理性，如解法所产生的误差能否满足精度要求：
解法是否稳定、是否收敛及熟练的速度等。
而且还能减少大量的人工计算。
由于工程实际中所遇到的数学模型求解过程迭代次数很多，计算量很大，所以需要借助如MATLAB，C++，VB，JAVA的辅助软件来解决，得到一个满足误差限的解。
本文所计算题目，均采用MATLAB进行编程，MATLAB被称为第四代计算机语言，利用其丰富的函数资源，使编程人员从繁琐的程序代码中解放出来
MATLAB最突出的特点就是简洁，它用更直观的、符合人们思维习惯的代码。
它具有以下优点：
1友好的工作平台和编程环境。
MATLAB界面精致，人机交互性强，操作简单。
2简单易用的程序语言。
MATLAB是一个高级的矩阵/阵列语言，包含控制语言、函数、数据结构，具有输入、输出和面向对象编程特点。
用户可以在命令窗口中将输入语句与执行命令同步，也可以先编好一个较大的复杂的应用程序（M文件）后再一起运行。
3强大的科学计算机数据处理能力。
包含大量计算算法的集合，拥有600多个工程中要用到的数学运算函数。
4出色的图像处理功能，可以方便地输出二维图像，便于我们绘制函数图像。
目 录
1.1实验目的4
1.2实验原理和方法4
1.3实验结果5
1.3.1最佳平方逼近法5
1.3.2拉格朗日插值法7
1.3.3对比8
2.2实验原理和方法10
2.3实验结果10
2.3.3第三问11
3.2实验原理和方法12
3.3实验结果12
1第一题
某过程涉及两变量x和y，拟分别用插值多项式和多项式拟合给出其对应规律的近似多项式，已知xi与yi之间的对应数据如下：
1
2
3
4
5
6
7
8
9
10
x
1
2
3
4
5
6
7
8
9
10
y
34.6588
40.3719
14.6448
-14.2721
-13.3570
24.8234
75.2795
103.5743
97.4847
78.2392
请用次数分别为3，4，5，6的多项式拟合并给出最好近似结果f（x）。
请用插值多项式给出最好近似结果。
1.1实验目的：
学习逼近和插值的原理和编程方法，由给出的已知点构造多项式，在某个范围内近似代替已知点所代表的函数，以便于简化对未知函数的各种计算。
1.2试验原理和方法：
实验原理：
拉格朗日插值法中先构造插值基础函数：
，然后构造出拉格朗日多项式：
。
最佳平方逼近中，设逼近函数
，逼近函数和真实函数之差
，
，即：
，根据最小二乘准则令
，可以得到
。
实验方法：
逼近法采用最佳平方逼近，依据最小二乘原则：
，由已知条件采用离散型。
插值法采用拉格朗日插值法。
在逼近法中，由于是离散型的，所以法方程系数阵设计成求和。
分别求出3、4、5、6次的多项式，逼近结果和真实值有一定差距，最小二乘正是让这些差距达到最小，理论上多项式次数越高结果和真实值差距越小。
拉格朗日插值法中“la=la*（p-x（j））/（x（k）-x（j））”语句实现的是我们通常书写的连乘形式拉格朗日插值多项式，但是表示不方便，而如果用“s=collect（s）”函数将其展开成降幂排列多项式以后，由于余项问题结果会和原本的多项式有偏差，这种偏差随着x的增大而增大。
求出多项式后和题目中给出的参考点进行比较。
最后，选择六次最佳平方逼近多项式和拉格朗日插值多项式（九次）进行比较，选取xi=a+ih=1+0.2*i（i=0,1,
45），分别绘制两者的图像进行比较。
1.3试验结果
1.3.1最佳平方逼近法
三次多项式：
-1.033*x^3+19.33*x^2-94.48*x+131.8
拟合结果：
1
2
3
4
5
6
7
8
9
10
x
1
2
3
4
5
6
7
8
9
10
y
55.6170
11.8960
-5.5610
-2.9520
13.5250
37.6720
63.2910
84.1840
94.1530
87.0000
四次多项式：
-0.3818*x^4+7.368*x^3-42.14*x^2+73.53*x+0.745
拟合结果：
1
2
3
4
5
6
7
8
9
10
x
1
2
3
4
5
6
7
8
9
10
y
39.1212
32.0802
10.0852
-5.5638
-2.7300
21.5602
61.1172
100.5882
115.4572
72.0450
五次多项式：
0.09807*t^5-3.079*t^4+34.5*t^3-163.5*t^2+304.7*t-139.5
拟合结果：
1
2
3
4
5
6
7
8
9
10
x
1
2
3
4
5
6
7
8
9
10
y
33.2191
45.7742
9.0320
-16.5003
-8.9063
26.9083
70.9835
100.0738
99.4164
74.5000
六次多项式：
0.01936*t^6-0.5408*t^5+5.114*t^4-16.9*t^3-0.867*t^2+66.38*t-18.7
拟合结果：
1
2
3
4
5
6
7
8
9
10
x
1
2
3
4
5
6
7
8
9
10
y
34.5056
41.1494
13.2700
-13.9486
-12.2250
23.7114
73.9500
105.1694
96.3456
78.4000
对比可知，六次多项式拟合结果最好。
1.3.2拉格朗日插值法
插值多项式5.353*10^（-5）*x^9-0.003088*x^8+0.07229*x^7-0.8792*x^6+5.932*x^5-22.41*x^4+50.11*x^3-86.47*x^2+113.5*x-25.2
注：
此多项式为拉格朗日多项式的近似式，当x=10的时候偏差可以达到23以上。
对比数据：
1
2
3
4
5
6
7
8
9
x
1.5000
1.9000
2.3000
2.7000
3.1000
3.5000
3.9000
4.3000
4.7000
y
42.1498
41.4620
35.1182
24.3852
11.2732
-1.7813
-12.3006
-18.1566
-17.9069
10
11
12
13
14
15
16
17
x
5.1000
5.5000
5.9000
6.3000
6.7000
7.1000
7.5000
7.9000
y
-11.0226
2.0284
19.8549
40.3626
61.0840
79.5688
93.7700
102.3677
插值结果：
1
2
3
4
5
6
7
8
9
x
1.5000
1.9000
2.3000
2.7000
3.1000
3.5000
3.9000
4.3000
4.7000
y
42.3840
41.4947
35.0742
24.3601
11.2792
-1.7683
-12.2977
-18.1626
-17.9118
10
11
12
13
14
15
16
17
x
5.1000
5.5000
5.9000
6.3000
6.7000
7.1000
7.5000
7.9000
y
-11.0210
2.0333
19.8565
40.3584
61.0794
79.5709
93.7788
102.3713
其中红点表示参考点。
1.3.3比较
选取xi=a+ih=1+0.2*i（i=0,1,
45），分别绘制六次多项式拟合和拉格朗日插值结果图：
其中绿线表示拉格朗日插值多项式图像，蓝线表示六次多项式拟合图像。
两者效果近似但后者比前者低三次。
2第二题
用雅格比法与高斯－赛德尔迭代法解下列方程组Ax=b1或Ax=b2，研究其收敛性。
上机验证理论分析是否正确，比较它们的收敛速度，观察右端项对迭代收敛有无影响。
（1）A行分别为A1=[6,2,-1],A2=[1,4,-2],A3=[-3,1,4]；b1=[-3,2,4]T；b2=[100,-200,345]T。
（2）A行分别为A1=[1,0,8,0.8],A2=[0.8,1,0.8],A3=[0.8,0.8,1]；b1=[3,2,1]T；b2=[5,0,-10]T。
（3）A行分别为A1=[1,3],A2=[-7,1]；b1=[4,6]T。
2.1试验目的
学习jacobi迭代法和GuassSeidel迭代法的原理和编程方法，研究方程组系数阵和右边项对方程的解及其收敛性的影响，判断迭代法的收敛条件。
2.2实验原理和方法
实验原理：
将方程组系数阵A分解为
其中D为对角阵，L为减去D的下三角阵，U为减去D的上三角阵。
Jacobi迭代法中构造如下迭代公式：
而Gauss-Seidel迭代法的迭代公式为：
初始值直接选取为0。
在判断其收敛性时，分别求解其迭代矩阵的谱半径
，
，
为迭代矩阵的特征值。
实验方法：
分别编写jacobi迭代及其收敛判别函数和Seidel迭代及其收敛判别函数。
如果在初试迭代步数之内还未收敛就进行收敛判别，收敛判别的依据是迭代矩阵的谱半径是否小于1。
比较同一方程组的jacobi迭代法和Seidel迭代法的结果是否相同，在达到精度要求后比较两种方法的迭代次数，比较哪一个的效率更高。
比较方程组系数阵和等号右边的变化会对方程的解和收敛速度造成什么影响。
如果迭代不收敛，那么考虑为什么不收敛，如果把方程组系数阵进行强对角占优处理，是否会收敛。
2.3实验结果
规定误差界：
1e-4
2.3.1第一问
.
由jacobi迭代法求得
，设定迭代20次，实际迭代16次，精度为9.4022e-005。
由seidel迭代法求得
，设定迭代20次，实际迭代10次，精度为：
9.0769e-005
，
由jacobi迭代法求得
，设定迭代40次，实际迭代23次，精度为6.9948e-005。
由Seidel迭代法求得
，设定迭代20次，实际迭代15次，精度为8.6384e-005
通过对比可知：
1、Seidel迭代的收敛速度明显高于jacobi迭代。
2、b矩阵对收敛速度和误差精度有影响，b中元素较大时会放慢收敛速度并加大误差。
2.3.2第二问
，
由jacobi迭代法求解，100次迭代尚且不能达到精度。
此时调用jacobi迭代法的收敛判别函数，求得特征值为：
，
，迭代不收敛。
由于Seidel迭代法求解
，迭代次数31，精度为8.7826e-005
，
Jacobi迭代不收敛。
Seldel迭代法求得
，迭代次数38，精度为8.4552e-005。
比较
得知A矩阵元素如果相差很小，迭代次数会大幅增加，综合比较
可知b矩阵元素如果相差很大会增加迭代次数。
2.3.3第三问试验结果和讨论
，
此时
，
，jacobi迭代法和Seidel迭代法都不收敛。
如果交换A中行的顺序，得到
，用jacobi迭代计算，迭代8次，解得
。
用Seidel迭代法计算，只需迭代5次，得到
，精度为2.6326e-005。
此时
，从此可以看出收敛速度的快慢。
3第三题
给定函数
，及节点
，求其三次样条插值多项式（可取I型或II型边界条件），并画图及与
的图形进行比较分析。
注：
涉及到线性方程组求解问题需采用适当的求解算法。
3.1实验目的
学习三弯矩法的原理和编程方法，对比原函数和三次样条插值的结果。
3.2实验原理和方法
实验原理：
记
，
，
，
，给出插值两端处的二阶导数
。
组成递推式：
由于系数阵按行对角占优，方程组存在唯一确定解，可以使用高斯列主元消去法来解方程。
最后将各个参数带入样条函数
即可求得样条函数。
实验方法
由于在本题中x（i+1）-x（i）=1,所以h（i）=1。
在编程中直接将h设置成常数，简化了运算。
首先求解μ、λ、g，然后列出方阵求解M（i），在求解方程组的过程中采用列主元素高斯消去法，分为消元和回代两个过程，编写这两个函数，解出除了两端的M，而两端点的M值等于两端点的函数二阶导数值。
编写函数求出样条函数的系数，然后求出方程，对于三弯矩法三次样条函数，如果有n个点，则有n-1个样条函数，除了两端需要求解n-2个M值，即解n-2阶方程组。
在表达样条函数的时候采用if语句，对不同的区间进行划分，然后细分[-5,5]这个区间，间隔0.1将其分为100份，这样可以体现出连续性，此时绘图对比三次样条函数和原函数。
本题中原函数的二阶导不是计算机解出的没有编写相关程序。
本题中求解的样条函数，MATLAB系统自动的将公因式提取，并且合并同类项，所以表达出的函数并不整齐规律，为了更好地体现三次样条函数的结构和性质，我专门手写了规整的样条函数。
3.3实验结果
三弯矩法求解
代入x=[-5-4-3-2-1012345]
y=[0.03850.05880.10000.20000.50001.00000.50000.20000.10000.05880.0385]
M=[0.01410.05990.09930.7431-1.87150.74310.09930.05990.0141]
求得样条函数系数阵为：
S=
0.00140.00240.03710.0565
0.00240.01000.05650.0900
0.01000.01650.09000.1835
0.01650.12380.18350.3762
0.1238-0.31190.37621.3119
-0.31190.12381.31190.3762
0.12380.01650.37620.1835
0.01650.01000.18350.0900
0.01000.00240.09000.0565
0.00240.00140.05650.0371
样条函数：
（97*X）/5000-（7*（X+4）^3）/5000+（3*（X+5）^3）/1250+1341/10000
（67*X）/2000-（3*（X+3）^3）/1250+（X+4）^3/100+381/2000
（187*X）/2000-（X+2）^3/100+（33*（X+3）^3）/2000+741/2000
（1927*X）/10000-（33*（X+1）^3）/2000+（619*（X+2）^3）/5000+5689/10000
（9357*X）/10000-（3119*（X+1）^3）/10000-（619*X^3）/5000+13119/10000
（3199*（X-1）^3）/10000-（9357*X）/10000+（619*X^3）/5000+13119/10000
（33*（X-1）^3）/2000-（1927*X）/10000-（619*（X-2）^3）/5000+5689/10000
（X-2）^3/100-（187*X）/2000-（33*（X-3）^3）/2000+741/2000
（3*（X-3）^3）/1250-（67*X）/2000-（X-4）^3/100+381/2000
（7*（X-4）^3）/5000-（97*X）/5000-（3*（X-5）^3）/1250+1341/10000
写出标准形式的样条函数：
S1（x）=0.0014*（-4-x）^3+0.0024*（x+5）^3+0.0371*（-4-x）+0.0565*（x+5）x[-5,-4]
S2（x）=0.0024*（-3-x）^3+0.0100*（x+4）^3+0.0565*（-3-x）+0.0900*（x+4）x[-4,-3]
S3（x）=0.0100*（-2-x）^3+0.0165*（x+3）^3+0.0900*（-2-x）+0.1835*（x+3）x[-3,-2]
S4（x）=0.0165*（-1-x）^3+0.1238*（x+2）^3+0.1835*（-1-x）+0.3762*（x+2）x[-2,-1]
S5（x）=0.1238*（-x）^3-0.3119*（x+1）^3+0.3762*（-x）+1.3119（x+1）x[-1,0]
S6（x）=-0.3119*（1-x）^3+0.1238*x^3+1.3119*（1-x）+0.3762*xx[0,1]
S7（x）=0.1238*（2-x）^3+0.0165*（x-1）^3+0.3762*（2-x）+0.1835*（x-1）x[1,2]
S8（x）=0.0165*（3-x）^3+0.0100*（x-2）^3+0.1835*（3-x）+0.0900*（x-2）x[2,3]
S9（x）=0.0100*（4-x）^3+0.0024*（x-3）^3+0.0900*（4-x）+0.0565*（x-3）x[3,4]
S10（x）=0.0024*（5-x）^3+0.0014*（x-4）^3+0.0565*（5-x）+0.0371*（x-4）x[4,5]
输入：
a=-5:
0.1:
5;
fori=1:
length（a）
b（i）=f（a（i））;
end
求得原函数和样条函数的对比图：
4MATLAB程序
第一题：
离散型最佳平方逼近函数
functionf=squar_approx_ls（x,y,n）
symst;
N=length（x）;
P=zeros（n+1）;
Q=zeros（n+1,1）;
a=0;
fori=0:
n
forj=0:
n
fork=1:
N
a=a+x（k）^（i+j）;
end
P（i+1,j+1）=a;
a=0;
end
end
b=0;
fori=0:
n
fork=1:
N
b=b+y（k）*x（k）^i
end
Q（i+1,1）=b;
b=0;
end
s=inv（P）*Q;
f=0;
fori=1:
n+1
f=f+s（i）*t^（i-1）;
simplify（f）;
end
f=collect（f）;
f=vpa（f,4）;
拉格朗日插值函数
functions=Lagrange（x,y,x0）
symsp;
n=length（x）;
s=0;
for（k=1:
n）
la=y（k）;
for（j=1:
k-1）
la=la*（p-x（j））/（x（k）-x（j））;
end;
for（j=k+1:
n）
la=la*（p-x（j））/（x（k）-x（j））;
end;
s=s+la;
simplify（s）;
end
if（nargin==2）
s=collect（s）;
s=vpa（s,4）;
else
m=length（x0）;
fori=1:
m
temp（i）=subs（s,'p',x0（i））;
end
s=temp;
end
第二题
Jacobi迭代法函数
functionx=jacobi（A,b,P,delta,n）
N=length（b）;
fork=1:
n
forj=1:
N
x（j）=（b（j）-A（j,[1:
j-1,j+1:
N]）*P（[1:
j-1,j+1:
N]））/A（j,j）;
end
err=abs（norm（x'-P））;
P=x';
if（errbreak;
end
P
end
x=x';
k,err
jacobi迭代收敛判别函数：
functions=liansanxing（A）
[N,N]=size（A）;
D=zeros（N）;
LU=zeros（N）;
fori=1:
N
D（i,i）=A（i,i）;
end
LU=A-D;
p=-inv（D）*LU;
[V,s]=eig（p）;
Seidel迭代法函数
functionx=Seidel（A,b,P,delta,n）
N=length（b）;
fork=1:
n
forj=1:
N
ifj==1
x
（1）=（b
（1）-A（1,2:
N）*P（2:
N））/A（1,1）;
elseifj==N
x（N）=（b（N）-A（N,1:
N-1）*（x（1:
N-1））'）/A（N,N）;
else
x（j）=（b（j）-A（j,1:
j-1）*x（1:
j-1）-A（j,j+1:
N）*P（j+1:
N））/A（j,j）;
end
end
err=abs（norm（x'-P））;
P=x';
if（errbreak;
end
end
P
x=x';
k,err
seidel迭代收敛判别函数：
functions=liansanxing2（A）
[N,N]=size（A）;
D=zeros（N）;
L=zeros（N）;
U=zeros（N）;
fori=2:
N
forj=1:
i-1
L（i,j）=A（i,j）;
end
end
fori=1:
N
D（i,i）=A（i,i）;
end
U=A-D-L;
p=-inv（D+L）*U;
[V,s]=eig（p）;
第三题
求解Mi函数
functions=sanwanju（x,y,y0d,ynd）
n=length（x）;
h=zeros（n-1,1）;
u=zeros（n-2,1
升级会员 