数值分析MATLAB上机实验
数值分析实习报告
姓名：gestepoA
学号：201*******
班级：***班
序言
随着计算机技术的迅速发展，数值分析在工程技术领域中的应用越来越广泛，并且成为数学与计算机之间的桥梁。要解决工程问题，往往需要处理很多数学模型，不仅要研究各种数学问题的数值解法，同时也要分析所用的数值解法在理论上的合理性，如解法所产生的误差能否满足精度要求：解法是否稳定、是否收敛及熟练的速度等。而且还能减少大量的人工计算。
由于工程实际中所遇到的数学模型求解过程迭代次数很多，计算量很大，所以需要借助如MATLAB，C++，VB，JAVA的辅助软件来解决，得到一个满足误差限的解。本文所计算题目，均采用MATLAB进行编程，MATLAB被称为第四代计算机语言，利用其丰富的函数资源，使编程人员从繁琐的程序代码中解放出来MATLAB最突出的特点就是简洁，它用更直观的、符合人们思维习惯的代码。它具有以下优点：
1友好的工作平台和编程环境。MATLAB界面精致，人机交互性强，操作简单。
2简单易用的程序语言。MATLAB是一个高级的矩阵/阵列语言，包含控制语言、函数、数据结构，具有输入、输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编好一个较大的复杂的应用程序（M 文件）后再一起运行。
3强大的科学计算机数据处理能力。包含大量计算算法的集合，拥有600多个工程中要用到的数学运算函数。
4出色的图像处理功能，可以方便地输出二维图像，便于我们绘制函数图像。
目录
1 第一题 (4)
1.1 实验目的 (4)
1.2 实验原理和方法 (4)
1.3 实验结果 (5)
1.3.1 最佳平方逼近法 (5)
1.3.2 拉格朗日插值法 (7)
1.3.3 对比 (8)
2 第二题 (9)
2.1实验目的 (9)
2.2 实验原理和方法 (10)
2.3 实验结果 (10)
2.3.1 第一问 (10)
2.3.2 第二问 (11)
2.3.3 第三问 (11)
3 第三题 (12)
3.1实验目的 (12)
3.2 实验原理和方法 (12)
3.3 实验结果 (12)
4 MATLAB程序 (14)
1第一题
某过程涉及两变量x 和y ，拟分别用插值多项式和多项式拟合给出其对应规律的
588
719
448 721
570
234
795 743
847
392
⑴请用次数分别为3，4，5，6的多项式拟合并给出最好近似结果f(x)。 ⑵请用插值多项式给出最好近似结果。
1.1实验目的：
学习逼近和插值的原理和编程方法，由给出的已知点构造多项式，在某个范围内近似代替已知点所代表的函数，以便于简化对未知函数的各种计算。
1.2试验原理和方法：
实验原理：
拉格朗日插值法中先构造插值基础函数：l k (x )=∏x?x i
x
k ?x i
n j=0j≠k
(k =0,1,2,?,n )，然后构
造出拉格朗日多项式：p n (x )=∑(∏x?x i
x
k ?x i
n j=0j≠k
)n k=0f (x k )。
最佳平方逼近中，设逼近函数P n (x )=a 0+a 1x +a 2x 2+?+a n x n ，逼近函数和真实函数之差r =P n (x )?y ，[r 1r 2?r n ]=[11?
x 1
x 2
????x 1n x 2n ?
1x n
?x n
n ][a 0a
1?a n ]?[y 1
y 2?y n ]，即：r =Xa ?Y ，根据最小二乘准则令∑r i 2n i=0=min ，可以得到a =(X T X )?1X T
Y 。
实验方法：
逼近法采用最佳平方逼近，依据最小二乘原则：∑r i 2
n i=0=min ，由已知条件采用离散型。插值法采用拉格朗日插值法。
在逼近法中，由于是离散型的，所以法方程系数阵设计成求和。分别求出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
拟合结果：
0000000000
四次多项式：- 0.3818*x^4 + 7.368*x^3 - 42.14*x^2 + 73.53*x+ 0.745
拟合结果：
y 39.121
232.080
2
10.085
2
-5.563
8
-2.730
21.560
2
61.117
2
100.588
2
115.457
2
72.045
五次多项式：0.09807*t^5 - 3.079*t^4 + 34.5*t^3 - 163.5*t^2 + 304.7*t - 139.5
拟合结果：
12345678910 x 1 2 3 4 5 6 7 8 9 10
y 33.219
120-16.500
3
-8.906
335840
六次多项式：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拟合结果：
y 34.505
640-13.94
86
-12.22
50409460
对比可知，六次多项式拟合结果最好。
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以上。
对比数据：
插值结果：
054493
其中红点表示参考点。
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分解为A=D+L+U,其中D为对角阵，L为减去D的下三角阵，U 为减去D的上三角阵。
Jacobi迭代法中构造如下迭代公式：
x(k+1)=?D?1(L+U)x(k)+D?1b
而Gauss-Seidel迭代法的迭代公式为：
x(k+1)=?(D+L)?1Ux(k)+(D+L)?1b
初始值直接选取为0。
在判断其收敛性时，分别求解其迭代矩阵的谱半径ρ(G)，ρ(G)=max
1≤i≤≤n
|λi|，λi为迭代矩阵的特征值。
实验方法：
分别编写jacobi迭代及其收敛判别函数和Seidel迭代及其收敛判别函数。如果在初试迭代步数之内还未收敛就进行收敛判别，收敛判别的依据是迭代矩阵的谱半径是否小于1。
比较同一方程组的jacobi迭代法和Seidel迭代法的结果是否相同，在达到精度要求后比较两种方法的迭代次数，比较哪一个的效率更高。
比较方程组系数阵和等号右边的变化会对方程的解和收敛速度造成什么影响。
如果迭代不收敛，那么考虑为什么不收敛，如果把方程组系数阵进行强对角占优处理，是否会收敛。
2.3实验结果
规定误差界：1e-4
2.3.1第一问
①A=[
62?1
1
?3
4
1
?2
4
],b=[
?3
2
4
].
由jacobi迭代法求得x=[ ?0.7273
0.8081
0.2525
]，设定迭代20次，实际迭代16次，精度为9.4022e-005。
由seidel迭代法求得x=[ ?0.7272
0.8081
0.2525
]，设定迭代20次，实际迭代10次，精度为：
9.0769e-005
②A=[62?1 1
?34
1
?2
4
]，b=[
100
?200
345
]
由jacobi 迭代法求得x =[ 36.3636
?2.0707114.0404]，设定迭代40次，实际迭代23次，精度为
6.9948e -005。
由Seidel 迭代法求得x =[ 36.3637
?2.0707114.0404]，设定迭代20次，实际迭代15次，精度为
8.6384e -005
③通过对比可知：1、Seidel 迭代的收敛速度明显高于jacobi 迭代。2、b 矩阵对收敛速度和误差精度有影响，b 中元素较大时会放慢收敛速度并加大误差。 2.3.2第二问
① A =[10.80.80.8
0.810.80.81
]，b =[3
21] 由jacobi 迭代法求解，100次迭代尚且不能达到精度。此时调用jacobi 迭代法的收敛判别函数，求得特征值为：λ1=?1.6，λ2、3=0.8，ρ(G J )=1.6>1，迭代不收敛。
由于Seidel 迭代法求解x =[ 5.7691
0.7693 ?4.2307
]，迭代次数31，精度为8.7826e-005
② A =[10.80.80.8
0.810.80.81
]，b =[5
0?10] Jacobi 迭代不收敛。
Seldel 迭代法求得x =[ 32.6922
7.6922 ?42.3076
] ，迭代次数38，精度为8.4552e -005。
比较①②得知A 矩阵元素如果相差很小，迭代次数会大幅增加，综合比较⑴⑵可知b 矩阵元素如果相差很大会增加迭代次数。 2.3.3第三问试验结果和讨论
A =[13
?71]，b =[46
]
此时ρ(G J )=4.5826>1，ρ(G G )=21>1，jacobi 迭代法和Seidel 迭代法都不收敛。
如果交换A 中行的顺序，得到[?7113]，用jacobi 迭代计算，迭代8次，解得x =[?0.6364
1.5454]。
用Seidel 迭代法计算，只需迭代5次，得到x =[ ?0.6364
1.5455]，精度为
2.6326e -005。此时ρ(G J )=
0.2182>ρ(G G )=0.0476，从此可以看出收敛速度的快慢。
3第三题
给定函数2
1
()5,15f x x
x -+≤≤=
，及节点50,1,10,i i x i ==-+，求其三次样条插值多项式（可取I 型或II 型边界条件），并画图及与()f x 的图形进行比较分析。注：涉及到线性方程组求解问题需采用适当的求解算法。
3.1实验目的
学习三弯矩法的原理和编程方法，对比原函数和三次样条插值的结果。
3.2实验原理和方法
实验原理：
记h i =x i ?x i?1，μi =
?i ?i +?i+1
，λi =1?μi ，g i =
6?i +?i+1
(
y i+1?y i ?i+1
?
y i ?y i?1
?i
)，给出插值
两端处的二阶导数S ′′(x 0)=y ′′(x 0)=M 0 S ′′(x n )=y ′′(x n )=M n 。组成递推式：
μi M i?1+2M i +λi M i+1=g i
由于系数阵按行对角占优，方程组存在唯一确定解，可以使用高斯列主元消去法来解方程。最后将各个参数带入样条函数
S i (x )=M i?1(x i ?x )36?i +M i (x ?x i?1)36?i +(y i?1?M i?16?i 2)x i ?x ?i +(y i ?M i 6?i 2)x ?x i?1
?i
即可求得样条函数。
实验方法
由于在本题中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 -1 0 1 2 3 4 5]
y=f(x)=
1 1+x2
f′′(?5)=f′′(5)≈0.00842
y=[0.0385 0.0588 0.1000 0.2000 0.5000 1.0000 0.5000 0.2000 0.1000 0.0588
0.0385]
M=[0.0141 0.0599 0.0993 0.7431 -1.8715 0.7431 0.0993 0.0599 0.0141]
求得样条函数系数阵为：
S=
0.0014 0.0024 0.0371 0.0565
0.0024 0.0100 0.0565 0.0900
0.0100 0.0165 0.0900 0.1835
0.0165 0.1238 0.1835 0.3762
0.1238 -0.3119 0.3762 1.3119
-0.3119 0.1238 1.3119 0.3762
0.1238 0.0165 0.3762 0.1835
0.0165 0.0100 0.1835 0.0900
0.0100 0.0024 0.0900 0.0565
0.0024 0.0014 0.0565 0.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*x x[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;
for i=1:length(a)
end
求得原函数和样条函数的对比图：
4 MATLAB程序
第一题：
离散型最佳平方逼近函数
function f=squar_approx_ls(x,y,n) syms t;
N=length(x);
P=zeros(n+1);
Q=zeros(n+1,1);
a=0;
for i=0:n
for j=0:n
for k=1:N
a=a+x(k)^(i+j);
end
P(i+1,j+1)=a;
a=0;
end
end
b=0;
for i=0:n
for k=1:N
b=b+y(k)*x(k)^i
Q(i+1,1)=b;
b=0;
end
s=inv(P)*Q;
f=0;
for i=1:n+1
f=f+s(i)*t^(i-1);
simplify(f);
end
f=collect(f);
f=vpa(f,4);
拉格朗日插值函数
function s=Lagrange(x,y,x0)
syms p;
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);
for i=1:m
temp(i)=subs(s,'p',x0(i));
end
s=temp;
end
第二题
Jacobi迭代法函数
function x=jacobi(A,b,P,delta,n) N=length(b);
for k=1:n
for j=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(err break; end P end x=x'; k,err jacobi迭代收敛判别函数： function s=liansanxing(A) [N,N]=size(A); D=zeros(N); LU=zeros(N); for i=1:N D(i,i)=A(i,i); end LU=A-D; p=-inv(D)*LU; [V,s]=eig(p); Seidel迭代法函数 function x=Seidel(A,b,P,delta,n) N=length(b); for k=1:n for j=1:N if j==1 x(1)=(b(1)-A(1,2:N)*P(2:N))/A(1,1); elseif j==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(err break; end end P x=x'; k,err seidel迭代收敛判别函数： function s=liansanxing2(A) [N,N]=size(A); D=zeros(N); L=zeros(N); U=zeros(N); for i=2:N for j=1:i-1 L(i,j)=A(i,j); end end for i=1:N D(i,i)=A(i,i); end U=A-D-L; p=-inv(D+L)*U; [V,s]=eig(p); 第三题 求解Mi函数 function s=sanwanju(x,y,y0d,ynd) n=length(x); h=zeros(n-1,1); u=zeros(n-2,1); v=zeros(n-2,1); g=zeros(n-2,1); for k=1:n-1 h(k)=x(k+1)-x(k); end for k=1:n-2 u(k)=h(k)/(h(k)+h(k+1)); v(k)=1-u(k); g(k)=(6/(h(k)+h(k+1)))*(((y(k+2)-y(k+1))/h(k+1))-((y(k+1)-y(k))/h (k))); end A=zeros(n-2); for i=1:n-2 A(i,i)=2; end for i=1:n-3 A(i,i+1)=v(i); end for i=2:n-2 A(i,i-1)=u(i); end B=zeros(n-2,1); B(1)=g(1)-u(1)*y0d; B(n-2)=g(n-2)-v(n-2)*ynd; for i=2:n-3 B(i)=g(i); end s=rowGauss(A,B); 高斯消去法函数 function x=rowGauss(A,b) [N,N]=size(A); x=zeros(1,N+1); Ag=[A,b]; for p=1:N-1 [y,j]=max(abs(Ag(p:N,p))); c=Ag(p,:); Ag(p,:)=Ag(j+p-1,:); Ag(j+p-1,:)=c; if Ag(p,p)==0 disp('奇异'); break; end for k=p+1:N m=Ag(k,p)/Ag(p,p); Ag(k,p:N+1)=Ag(k,p:N+1)-m*Ag(p,p:N+1); end end x=backsub(Ag(1:N,1:N),Ag(1:N,N+1)); 回代法函数 function x=backsub(A,b) n=length(b); x=zeros(n,1); x(n)=b(n)/A(n,n); for k=n-1:-1:1 x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k); end 基于本题的样条函数求解函数 function s=sanwanjufc(y,m,m0,mn) n=length(y); a=zeros(n-1,1); b=zeros(n-1,1); c=zeros(n-1,1); d=zeros(n-1,1); a(1)=m0/6; for i=2:n-1 a(i)=m(i-1)/6; end b(n-1)=mn/6; for i=1:n-2 b(i)=m(i)/6; end c(1)=y(1)-m0/6; for i=2:n-1 c(i)=y(i)-m(i-1)/6; end d(n-1)=y(n)-mn/6; for i=1:n-2 d(i)=y(i+1)-m(i)/6; end S=[a b c d]; S; for i=1:n-1 s(i)=a(i)*(x(i+1)-X)^3+b(i)*(X-x(i))^3+c(i)*(x(i+1)-X)+d(i)*(X-x( i)); simplify(s(i)); disp(s(i)); end 样条函数表达函数 function y=fd(x) if x>=-5&&x<=-4 y=0.0014*(-4-x)^3+0.0024*(x+5)^3+0.0371*(-4-x)+0.0565*(x+5); elseif x>=-4&&x<=-3 y=0.0024*(-3-x)^3+0.0100*(x+4)^3+0.0565*(-3-x)+0.0900*(x+4); elseif x>=-3&&x<=-2 y=0.0100*(-2-x)^3+0.0165*(x+3)^3+0.0900*(-2-x)+0.1835*(x+3); elseif x>=-2&&x<=-1 y=0.0165*(-1-x)^3+0.1238*(x+2)^3+0.1835*(-1-x)+0.3762*(x+2); elseif x>=-1&&x<=0 y=0.1238*(-x)^3-0.3119*(x+1)^3+0.3762*(-x)+1.3119*(x+1); elseif x>=0&&x<=1 y=-0.3119*(1-x)^3+0.1238*x^3+1.3119*(1-x)+0.3762*x elseif x>=1&&x<=2 y=0.1238*(2-x)^3+0.0165*(x-1)^3+0.3762*(2-x)+0.1835*(x-1); elseif x>=2&&x<=3 y=0.0165*(3-x)^3+0.0100*(x-2)^3+0.1835*(3-x)+0.0900*(x-2); elseif x>=3&&x<=4 y=0.0100*(4-x)^3+0.0024*(x-3)^3+0.0900*(4-x)+0.0565*(x-3); elseif x>=4&&x<=5 y=0.0024*(5-x)^3+0.0014*(x-4)^3+0.0565*(5-x)+0.0371*(x-4); end 数值分析实习报告 姓名：gestepoA 学号：201******* 班级：***班 序言 随着计算机技术的迅速发展，数值分析在工程技术领域中的应用越来越广泛，并且成为数学与计算机之间的桥梁。要解决工程问题，往往需要处理很多数学模型，不仅要研究各种数学问题的数值解法，同时也要分析所用的数值解法在理论上的合理性，如解法所产生的误差能否满足精度要求：解法是否稳定、是否收敛及熟练的速度等。而且还能减少大量的人工计算。 由于工程实际中所遇到的数学模型求解过程迭代次数很多，计算量很大，所以需要借助如MATLAB，C++，VB，JAVA的辅助软件来解决，得到一个满足误差限的解。本文所计算题目，均采用MATLAB进行编程，MATLAB被称为第四代计算机语言，利用其丰富的函数资源，使编程人员从繁琐的程序代码中解放出来MATLAB最突出的特点就是简洁，它用更直观的、符合人们思维习惯的代码。它具有以下优点： 1友好的工作平台和编程环境。MATLAB界面精致，人机交互性强，操作简单。 2简单易用的程序语言。MATLAB是一个高级的矩阵/阵列语言，包含控制语言、函数、数据结构，具有输入、输出和面向对象编程特点。用户可以在命令窗口中将输入语句与执行命令同步，也可以先编好一个较大的复杂的应用程序（M 文件）后再一起运行。 3强大的科学计算机数据处理能力。包含大量计算算法的集合，拥有600多个工程中要用到的数学运算函数。 4出色的图像处理功能，可以方便地输出二维图像，便于我们绘制函数图像。 目录 1 第一题 (4) 1.1 实验目的 (4) 1.2 实验原理和方法 (4) 1.3 实验结果 (5) 1.3.1 最佳平方逼近法 (5) 1.3.2 拉格朗日插值法 (7) 1.3.3 对比 (8) 2 第二题 (9) 2.1实验目的 (9) 2.2 实验原理和方法 (10) 2.3 实验结果 (10) 2.3.1 第一问 (10) 2.3.2 第二问 (11) 2.3.3 第三问 (11) 3 第三题 (12) 3.1实验目的 (12) 3.2 实验原理和方法 (12) 3.3 实验结果 (12) 4 MATLAB程序 (14) MATLAB数值计算 MATLAB数值计算 (1) 1创建矩阵 (3) 1.1直接输入 (3) 1.2向量 (3) 1.2.1linspace：线性分布 (3) 1.2.2冒号法 (3) 1.3函数创建 (4) 1.3.1eye:单位矩阵 (4) 1.3.2rand：随机矩阵 (4) 1.3.3zeros：全0矩阵 (4) 1.3.4ones：全1矩阵 (5) 2矩阵运算 (5) 2.1加减 (5) 2.1.1[M×N]±[M×N] (5) 2.2乘 (6) 2.2.1[M×N]*a (6) 2.2.2[M×N]*[N×M] (6) 2.3乘方 (7) 2.3.1[M×M]^a (7) 2.3.2a^[M×M] (7) 2.4特殊运算 (8) 2.4.1求逆inv (8) 2.4.2行列式det (8) 2.4.3特征值eig (8) 2.4.4转置'和.' (9) 2.4.5变形reshape (10) 2.4.6翻转rot90,fliplr,flipud (11) 2.4.7抽取diag,tril,triu (12) 2.5数组运算 (12) 2.5.1乘 (12) [M×N].*[M×N] (12) 2.5.2除 (13) [M×N]./[M×N] (14) [M×N].\[M×N] (14) 2.5.3乘方 (14) [M×N].^[M×N] (15) a.^[M×N] (15) 2.6除法 (15) 2.6.1求解线性方程组 (15) 3多项式 (16) 3.1系数表示法poly (16) 3.2求根roots (16) 3.3乘法conv (16) 3.4除法deconv (17) 3.5求值polyval (17) 3.6微分polyder (18) 电子科技大学电子工程学院标准实验报告（实验）课程名称MATLAB与数值分析 学生姓名：李培睿 学号：2013020904026 指导教师：程建 一、实验名称 《MATLAB与数值分析》第一次上机实验 二、实验目的 1. 熟练掌握矩阵的生成、加、减、乘、除、转置、行列式、逆、范数等运算 操作。（用.m文件和Matlab函数编写一个对给定矩阵进行运算操作的程序） 2. 熟练掌握算术符号操作和基本运算操作，包括矩阵合并、向量合并、符号 转换、展开符号表达式、符号因式分解、符号表达式的化简、代数方程的符号解析解、特征多项式、函数的反函数、函数计算器、微积分、常微分方程的符号解、符号函数的画图等。（用.m文件编写进行符号因式分解和函数求反的程序） 3. 掌握Matlab函数的编写规范。 4、掌握Matlab常用的绘图处理操作，包括：基本平面图、图形注释命令、 三维曲线和面的填充、三维等高线等。（用.m文件编写在一个图形窗口上绘制正弦和余弦函数的图形，并给出充分的图形注释） 5. 熟练操作MATLAB软件平台，能利用M文件完成MATLAB的程序设计。 三、实验内容 1. 编程实现以下数列的图像，用户能输入不同的初始值以及系数。并以x， y为坐标显示图像 x(n+1) = a*x(n)-b*(y(n)-x(n)^2); y(n+1) = b*x(n)+a*(y(n)-x(n)^2) 2. 编程实现奥运5环图，允许用户输入环的直径。 3. 实现对输入任意长度向量元素的冒泡排序的升序排列。不允许使用sort 函数。 四、实验数据及结果分析 题目一： ①在Editor窗口编写函数代码如下： Matlab作业3（数值分析） 机电工程学院（院、系）专业班组 学号姓名实验日期教师评定 1.计算多项式乘法(x2+2x+2)(x2+5x+4)。 答： 2. （1）将(x-6)(x-3)(x-8)展开为系数多项式的形式。（2）求解在x=8时多项 式(x-1)(x-2) (x-3)(x-4)的值。 答：（1） （2） 3. y=sin(x)，x从0到2π，?x=0.02π，求y的最大值、最小值、均值和标准差。 4.设ｘ＝[0.00.30.8 1.1 1.6 2.3]'，y=[0.500.82 1.14 1.25 1.35 1.40]'，试求二次多项式拟合系数，并据此计算x1=[0.9 1.2]时对应的y1。解：x=[0.0 0.3 0.8 1.1 1.6 2.3]'; %输入变量数据x y=[0.50 0.82 1.14 1.25 1.35 1.40]'; %输入变量数据y p=polyfit(x,y,2) %对x,y用二次多项式拟合,得到系数p x1=[0.9 1.2]; %输入点x1 y1=polyval(p,x1) %估计x1处对应的y1 p = -0.2387 0.9191 0.5318 y1 = a) 1.2909 5.实验数据处理：已知某压力传感器的测试数据如下表 p为压力值，u为电压值，试用多项式 d cp bp ap p u+ + + =2 3 ) ( 来拟 合其特性函数，求出a,b,c,d，并把拟合曲线和各个测试数据点画在同一幅图上。解： >> p=[0.0,1.1,2.1,2.8,4.2,5.0,6.1,6.9,8.1,9.0,9.9]; u=[10,11,13,14,17,18,22,24,29,34,39]; x=polyfit(p,u,3) %得多项式系数 t=linspace(0,10,100); y=polyval(x,t); %求多项式得值 plot(p,u,'*',t,y,'r') %画拟和曲线 x = 0.0195 -0.0412 1.4469 9.8267 第2章牛顿插值法实现 参考文献：[1]岑宝俊. 牛顿插值法在凸轮曲线修正设计中的应用[J]. 机械工程师,2009,10:54-55. 求牛顿插值多项式和差商的MA TLAB 主程序： function[A,C,L,wcgs,Cw]=newpoly(X,Y) n=length(X);A=zeros(n,n);A(:,1) =Y'; s=0.0;p=1.0;q=1.0;c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1)); end b=poly(X(j-1));q1=conv(q,b);c1=c1*j;q=q1; end C=A(n,n);b=poly(X(n));q1=conv(q1,b); for k=(n-1):-1:1 C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k); end L(k,:)=poly2sym(C);Q=poly2sym(q1); syms M wcgs=M*Q/c1;Cw=q1/c1; （1）保存名为newpoly.m 的M 文件 （2）输入MA TLAB 程序 >> X=[242,243,249,250]; >> Y=[13.681,13.526,13.098,13.095]; >> [A,C,L,wcgs,Cw]=newpoly(X,Y) 输出3阶牛顿插值多项式L 及其系数向量C 差商的矩阵A ，插值余项wcgs 及其 ) ()()1(ξ+n n f x R 的系数向量Cw 。 A = 13.6810 0 0 0 13.5260 -0.1550 0 0 13.0980 -0.0713 0.0120 0 13.0950 -0.0030 0.0098 -0.0003 C = 1.0e+003 * MATLAB 实验一 MATLAB 数值计算 试验报告说明： 1 做试验前请先预习，并独立完成试验和试验报告。 2 报告解答方式：将MATLAB 执行命令和最后运行结果从命令窗口拷贝到每题的题目下面，请将报告解答部分的底纹设置为灰色，以便于批阅。 3 在页眉上写清报告名称，学生姓名，学号，专业以及班级。 3 报告以Word 文档书写。 文档命名方式: 学号+姓名+_（下划线）+试验几.doc 如：110400220张三_试验1.doc 4 试验报告doc 文档以附件形式发送到maya_email@https://www.360docs.net/doc/082975118.html, 。凡文档命名不符合规范，或者发送方式不正确，不予登记。 5 每次试验报告的最后提交期限：下次试验课之前。 一 目的和要求 1 熟练掌握MATLAB 变量的使用 2 熟练掌握矩阵的创建 3 熟练掌握MATLAB 的矩阵和数组的运算 4 使用元胞数组和结构数组 二 试验内容 1 创建矩阵（必做） 1.1使用直接输入，from:step:to ，linspace ，logspace 等方式创建矩阵。 1.2 输入矩阵12342 46836 9 12a ?? ?= ? ?? ? 1.2-1）分别使用全下标和单下标达方式取出元素“8” >>a=[1 2 3 4;2 4 6 8;3 6 9 12] >> a(2,4) %全下标方式 >> a(11) % 单下标方式 1.2-2）分别用不同的方式从矩阵a 中取出子矩阵?? ??? ?1286 4 3 2 %方法一：全下标方式 a([2,3],[1 2 4]) %方法二：单下标方式 a([2 5 11;3 6 12]) % 方法三：利用逻辑向量 l1=logical([0 1 1]) 列主元法 function lianzhuyuan(A,b) n=input('请输入n:') %选择阶数A=zeros(n,n); %系数矩阵A b=zeros(n,1); %矩阵b X=zeros(n,1); %解X for i=1:n for j=1:n A(i,j)=(1/(i+j-1)); %生成hilbert矩阵A end b(i,1)=sum(A(i,:)); %生成矩阵b end for i=1:n-1 j=i; top=max(abs(A(i:n,j))); %列主元 k=j; while abs(A(k,j))~=top %列主元所在行 k=k+1; end for z=1:n %交换主元所在行a1=A(i,z); A(i,z)=A(k,z); A(k,z)=a1; end a2=b(i,1); b(i,1)=b(k,1); b(k,1)=a2; for s=i+1:n %消去算法开始m=A(s,j)/A(i,j); %化简为上三角矩阵 A(s,j)=0; for p=i+1:n A(s,p)=A(s,p)-m*A(i,p); end b(s,1)=b(s,1)-m*b(i,1); end end X(n,1)=b(n,1)/A(n,n); %回代开始 for i=n-1:-1:1 s=0; %初始化s for j=i+1:n s=s+A(i,j)*X(j,1); end X(i,1)=(b(i,1)-s)/A(i,i); end X 欧拉法 clc clear % 欧拉法 p=10; %贝塔的取值 T=10; %t取值的上限 y1=1; %y1的初值 r1=1; %y2的初值 %输入步长h的值 h=input('欧拉法please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0 break end S1=0:T/h; S2=0:T/h; S3=0:T/h; S4=0:T/h; i=1; % 迭代过程 for t=0:h:T Y=(exp(-t)); R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t); y=y1+h*(-y1); y1=y; r=r1+h*(y1-p*r1); r1=r; S1(i)=Y; S2(i)=R; S3(i)=y; S4(i)=r; i=i+1; end t=[0:h:T]; % 红线为解析解，'x'为数值解 plot(t,S1,'r',t,S3,'x') MATLAB 编程题库 1.下面的数据表近似地满足函数2 1cx b ax y ++=，请适当变换成为线性最小二乘问题，编程求最好的系数c b a ,,，并在同一个图上画出所有数据和函数图像. 625 .0718.0801.0823.0802.0687.0606.0356.0995 .0628.0544.0008.0213.0362.0586.0931.0i i y x ---- 解： x=[-0.931 -0.586 -0.362 -0.213 0.008 0.544 0.628 0.995]'; y=[0.356 0.606 0.687 0.802 0.823 0.801 0.718 0.625]'; A=[x ones(8,1) -x.^2.*y]; z=A\y; a=z(1); b=z(2); c=z(3); xh=-1:0.1:1; yh=(a.*xh+b)./(1+c.*xh.^2); plot(x,y,'r+',xh,yh,'b*') 2.若在Matlab工作目录下已经有如下两个函数文件，写一个割线法程序，求出这两个函数 10 的近似根，并写出调用方式： 精度为10 解： >> edit gexianfa.m function [x iter]=gexianfa(f,x0,x1,tol) iter=0; while(norm(x1-x0)>tol) iter=iter+1; x=x1-feval(f,x1).*(x1-x0)./(feval(f,x1)-feval(f,x0)); x0=x1;x1=x; end >> edit f.m function v=f(x) v=x.*log(x)-1; >> edit g.m function z=g(y) z=y.^5+y-1; >> [x1 iter1]=gexianfa('f',1,3,1e-10) x1 = 1.7632 iter1 = 6 >> [x2 iter2]=gexianfa('g',0,1,1e-10) x2 = 0.7549 iter2 = 8 实验一 误差分析 实验（病态问题） 实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出：考虑一个高次的代数多项式 )1.1() ()20()2)(1()(20 1∏=-=---=k k x x x x x p 显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动 )2.1(0 )(19=+x x p ε 其中ε是一个非常小的数。这相当于是对（）中19x 的系数作一个小的扰动。我们希望比较（）和（）根的差别，从而分析方程（）的解对扰动的敏感性。 实验内容：为了实现方便，我们先介绍两个Matlab 函数：“roots ”和“poly ”。 roots(a)u = 其中若变量a 存储n+1维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,,,+n a a a ，则输出u 的各分量是多项式方程 01121=+++++-n n n n a x a x a x a 的全部根；而函数 poly(v)b = 的输出b 是一个n+1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“poly ”是两个互逆的运算函数。 ;000000001.0=ess );21,1(zeros ve = ;)2(ess ve = ))20:1((ve poly roots + 上述简单的Matlab 程序便得到（）的全部根，程序中的“ess ”即是（）中的ε。 实验要求： （1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。 如果扰动项的系数ε很小，我们自然感觉（）和（）的解应当相差很小。计算中你有什么出乎意料的发现表明有些解关于如此的扰动敏感性如何 （2）将方程（）中的扰动项改成18x ε或其它形式，实验中又有怎样的现象 出现 （3）（选作部分）请从理论上分析产生这一问题的根源。注意我们可以将 方程（）写成展开的形式， ) 3.1(0 ),(1920=+-= x x x p αα 同时将方程的解x 看成是系数α的函数，考察方程的某个解关于α的扰动是否敏感，与研究它关于α的导数的大小有何关系为什么你发现了什么现象，哪些根关于α的变化更敏感 思考题一：（上述实验的改进） 在上述实验中我们会发现用roots 函数求解多项式方程的精度不高，为此你可以考虑用符号函数solve 来提高解的精确度，这需要用到将多项式转换为符号多项式的函数poly2sym,函数的具体使用方法可参考Matlab 的帮助。 Tutorial 6 实验报告 实验名称：Matlab数值计算 实验目的： 1、掌握数据统计与分析的方法； 2、掌握数据插值和曲线拟合的方法及其应用； 3、掌握多项式的常用运算。 实验内容： 1.利用randn函数生成符合正态分布的10×5随机矩阵A，进行如下操作： （1）求A的最大元素和最小元素； （2）求A的每行元素的和以及全部元素的和； （3）分别对A的每列元素按升序、每行元素按降序排列。 2.用3次多项式方法插值计算1-100之间整数的平方根。 3.某气象观测站测得某日6：00-18：00之间每隔2h的室内外温度（°C）如下表所示。 使用三次样条插值分别求出该日室内外6：30-17：30之间每隔2h各点的近似温度，并绘制插值后的温度曲线。 4.已知lgx在[1,101]区间10个整数采样点的函数值如下表所示， 试求lgx 的5次拟合多项式p(x)，并绘制lgx 和p(x)在[1,101]区间的函数曲线。 5. 有3个多项式(),(),()P x x x x P x x P x x x =+++=+=++4322 123 24522 3，试进行下列操作： （1） 求()()()()P x P x P x P x =+123。 （2） 求()P x 的根。 （3） 当x 取矩阵A 的每一元素时，求()P x 的值。其中： .....A --?? ??=?? ???? 112140752350525 6. 求函数在指定点的数值导数。 (),,f x x ==123 7. 用数值方法求定积分。 （1）I π =? 210的近似值。 （2）ln() x I dx x += +?1 22011 实验结果： 1. 1、%用牛顿法求f(x)=x-sin x 的零点，e=10^(-6) disp('牛顿法'); i=1; n0=180; p0=pi/3; tol=10^(-6); for i=1:n0 p=p0-(p0-sin(p0))/(1-cos(p0)); if abs(p-p0)<=10^(-6) disp('用牛顿法求得方程的根为') disp(p); disp('迭代次数为：') disp(i) break; end p0=p; end if i==n0&&~(abs(p-p0)<=10^(-6)) disp(n0) disp('次牛顿迭代后无法求出方程的解') end 2、disp('Steffensen加速'); p0=pi/3; for i=1:n0 p1=0.5*p0+0.5*cos(p0); p2=0.5*p1+0.5*cos(p1); p=p0-((p1-p0).^2)./(p2-2.*p1+p0); if abs(p-p0)<=10^(-6) disp('用Steffensen加速求得方程的根为') disp(p); disp('迭代次数为：') disp(i) break; end p0=p; end if i==n0&&~(abs(p-p0)<=10^(-6)) disp(n0) disp('次Steffensen加速后无法求出方程的解') end 1、%使用二分法找到方程 600 x^4 -550 x^3 +200 x^2 -20 x -1 =0 在区间[0.1,1]上的根， %误差限为 e=10^-4 disp('二分法') a=0.2;b=0.26; tol=0.0001; n0=10; fa=600*(a.^4)-550*(a.^3)+200*(a.^2)-20*a-1; for i=1:n0 p=(a+b)/2; fp=600*(p.^4)-550*(p.^3)+200*(p.^2)-20*p-1; if fp==0||(abs((b-a)/2) 实验6 Matlab数值计算 实验目的： 1、掌握数据统计与分析的方法； 2、掌握数据插值和曲线拟合的方法及其应用； 3、掌握多项式的常用运算。 实验内容： 1.利用randn函数生成符合正态分布的10×5随机矩阵A，进行如下操作： （1）求A的最大元素和最小元素； （2）求A的每行元素的和以及全部元素的和； （3）分别对A的每列元素按升序、每行元素按降序排列。 a = randn(10,5)+10; ma = max(max(a)) mi = min(min(a)) s = sum(a,2) sa = sum(sum(a)) p = sort(a) p1 = -sort(-a,2) 2.用3次多项式方法插值计算1-100之间整数的平方根。 f = sqrt(n); interp1(n,f,(1:100),'cubic') 3.某气象观测站测得某日6：00-18：00之间每隔2h的室内外温度（°C）如下表所示。 使用三次样条插值分别求出该日室内外6：30-17：30之间每隔2h 各点的近似温度，并绘制插值后的温度曲线。 n= 6:2:18; f1 = [18 20 22 25 30 28 24]; f2 = [15 19 24 28 34 32 30]; r = 6.5:2:17.5; w = interp1(n,f1,r,'spline'); w1 = interp1(n,f2,r,'spline'); subplot(211),plot(r,w) subplot(212),plot(r,w1) 4. 已知lgx 在[1,101]区间10个整数采样点的函数值如下表所示， 试求lgx 的5次拟合多项式p(x)，并绘制lgx 和p(x)在[1,101]区间的函数曲线。 x = linspace(1,101,10); y = log(x) /log(10); p = polyfit(x,y,5) y1 = polyval(p,x) plot(x,y,':o',x,y1,'-*') legend('sin(x)','fit') 5. 有3个多项式(),(),()P x x x x P x x P x x x =+++=+=++4 3 2 2 123245223，试进 行下列操作： （1） 求()()()()P x P x P x P x =+123。 （2） 求()P x 的根。 （3） 当x 取矩阵A 的每一元素时，求()P x 的值。其中： .....A --?? ? ?=?????? 11214075 2350 5 25 p1 = [1 2 4 0 5]; p2 = [0 0 0 1 2]; 第二讲MATLAB的数值分析 2-1矩阵运算与数组运算 矩阵运算和数组运算是MATLAB数值运算的两大类型，矩阵运算是按矩阵的运算规则进行的，而数组运算则是按数组元素逐一进行的。因此，在进行某些运算（如乘、除）时，矩阵运算和数组运算有着较大的差别。在MATLAB中，可以对矩阵进行数组运算，这时是把矩阵视为数组，运算按数组的运算规则。也可以对数组进行矩阵运算，这时是把数组视为矩阵，运算按矩阵的运算规则进行。 1、矩阵加减与数组加减 矩阵加减与数组加减运算效果一致，运算符也相同，可分为两种情况： （1）若参与运算的两矩阵（数组）的维数相同，则加减运算的结果是将两矩阵的对应元素进行加减，如 A=[1 1 1;2 2 2;3 3 3]; B=A; A+B ans= 2 2 2 4 4 4 6 6 6 （2）若参与运算的两矩阵之一为标量（1*1的矩阵），则加减运算的结果是将矩阵（数组）的每一元素与该标量逐一相加减，如 A=[1 1 1;2 2 2;3 3 3]; A+2 ans= 3 3 3 4 4 4 5 5 5 2、矩阵乘与数组乘 （1）矩阵乘 矩阵乘与数组乘有着较大差别，运算结果也完全不同。矩阵乘的运算符为“*”，运算是按矩阵的乘法规则进行，即参与乘运算的两矩阵的内维必须相同。设A、B为参与乘运算的 =A m×k B k×n。因此，参与运两矩阵，C为A和B的矩阵乘的结果，则它们必须满足关系C m ×n 算的两矩阵的顺序不能任意调换，因为A*B和B*A计算结果很可能是完全不一样的。如：A=[1 1 1;2 2 2;3 3 3]; B=A; A*B ans= 6 6 6 12 12 12 18 18 18 F=ones(1,3); G=ones(3,1); F*G ans 3 G*F ans= 1 1 1 1 1 1 1 1 1 （2）数组乘 数组乘的运算符为“.*”，运算符中的点号不能遗漏，也不能随意加空格符。参加数组乘运算的两数组的大小必须相等（即同维数组）。数组乘的结果是将两同维数组（矩阵）的对应元素逐一相乘，因此，A.*B和B.*A的计算结果是完全相同的，如： A=[1 1 1 1 1;2 2 2 2 2;3 3 3 3 3]; B=A; A.*B ans= 1 1 1 1 1 4 4 4 4 4 9 9 9 9 9 B.*A ans= 1 1 1 1 1 4 4 4 4 4 9 9 9 9 9 由于矩阵运算和数组运算的差异，能进行数组乘运算的两矩阵，不一定能进行矩阵乘运算。如 A=ones(1,3); B=A; A.*B ans= 1 1 1 A*A ???Error using= => 佛山科学技术学院 实 验 报 告 课程名称 数值分析 实验项目 数值积分 专业班级 机械工程 姓 名 余红杰 学 号 2111505010 指导教师 陈剑 成 绩 日 期 月 日 一、实验目的 1、 掌握程序的录入和matlab 的使用和操作； 2、 了解影响线性方程组解的精度的因素——方法与问题的性态。 3、 学会Matlab 提供的“\”的求解线性方程组。 二、实验要求 1、按照题目要求完成实验内容； 2、写出相应的Matlab 程序； 3、给出实验结果(可以用表格展示实验结果)； 4、分析和讨论实验结果并提出可能的优化实验。 5、写出实验报告。 三、实验步骤 1、用LU 分解及列主元高斯消去法解线性方程组 a)??????? ??=??????? ????????? ??----15900001.582012151526099999.2310 7104321x x x x ， 输出b Ax =中系数LU A =分解的矩阵L 和U ，解向量x 和)det(A ;用列主元法的行交换次序解向量x 和求)det(A ；比较两种方法所得结果。 2、用列主高斯消元法解线性方程组b Ax =。 （1）、???? ? ??=????? ??????? ??--11134.981.4987.023.116.427 .199.103.601.3321x x x （2）、???? ? ??=????? ??????? ??--11134.981.4990.023.116.427 .199.103.600.3321x x x 分别输出)det(,,A b A ，解向量x ，（1）中A 的条件数。分析比较（1）、（2）的计算结果 3、线性方程组b Ax =的A 和b 分别为 ??????? ??=1095791068565778710A ,?????? ? ??=31332332b 则解T x ),1,1,1,1(=. 用MATLAB 内部函数求)det(A 和A 的所有特征值和2)(A cond . 若令 ?????? ? ??=+98.99599.6989.998.585604.508.72.71.8710A A δ, 求解b x x A A =++))((δδ，输出向量x δ和2x δ，从理论结果和实际计算两方面分析线性 方程组b Ax =解的相对误差22/x x δ以及A 的相对误差 /A A δ的关系。 四、实验结果 1： %run311.m clc,clear; A = [10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2]; b = [8;5.90001;5;1]; %L U 分解 format short %小数点后四位，不然会受到后面的影响 [L U] = lu(A) %解方程组，输出A ，det(A) y = L\b; format long %小数点后15位显示 x = U\y 数值分析幂法与反幂法 matlab程序 随机产生一对称矩阵，对不同的原点位移和初值(至少取3个)分别使用幂法求计算矩阵的主特征值及主特征向量，用反幂法求计算矩阵的按模最小特征值及特征向量。 要求 1）比较不同的原点位移和初值说明收敛性 2）给出迭代结果，生成DOC文件。 3）程序清单，生成M文件。 解答： >> A=rand(5) %随机产生5*5矩阵求随机矩阵 A = 0.7094 0.1626 0.5853 0.6991 0.1493 0.7547 0.1190 0.2238 0.8909 0.2575 0.2760 0.4984 0.7513 0.9593 0.8407 0.6797 0.9597 0.2551 0.5472 0.2543 0.6551 0.3404 0.5060 0.1386 0.8143 >> B=A+A' %A矩阵和A的转置相加，得到随机对称矩阵B B = 1.4187 0.9173 0.8613 1.3788 0.8044 0.9173 0.2380 0.7222 1.8506 0.5979 0.8613 0.7222 1.5025 1.2144 1.3467 1.3788 1.8506 1.2144 1.0944 0.3929 0.8044 0.5979 1.3467 0.3929 1.6286 B=?? ????? ???? ?? ???6286.13929.03467.15979.08044 .03929.00944 .12144.18506 .13788.13467.12144.15025.17222.08613.05979.08506.17222.02380.09173.08044.03788.18613 .09173 .04187.1 编写幂法、反幂法程序： function [m,u,index,k]=pow(A,u,ep,it_max) % 求矩阵最大特征值的幂法，其中 % A 为矩阵； % ep 为精度要求，缺省为1e-5； % it_max 为最大迭代次数，缺省为100； % m 为绝对值最大的特征值； % u 为对应最大特征值的特征向量； % index ，当index=1时，迭代成功，当index=0时，迭代失败 if nargin<4 it_max=100; end if nargin<3 ep=1e-5; end n=length(A); index=0; k=0; m1=0; m0=0.01; % 修改移位参数，原点移位法加速收敛，为0时，即为幂法 I=eye(n) T=A-m0*I while k<=it_max v=T*u; [vmax,i]=max(abs(v)); m=v(i); u=v/m; if abs(m-m1) 基于MATLAB数值分析实验报告 班级：072115 姓名：李凯 学号：20111003943 实验二：矩阵与向量运算 实验目的：在MATLAB里，会对矩阵与向量进行加、减、数乘、求逆及矩阵特征值运算，以及矩阵的LU分解。 设A是一个n×n方阵，X是一个n维向量，乘积Y=AX可以看作是n维空间变换。如果能够找到一个标量λ，使得存在一个非零向量X，满足：AX=λX （3.1）则可以认为线性变换T（X）=AX将X映射为λX，此时，称X 是对应于特征值λ的特征向量。改写式（3.1）可以得到线性方程组的标准形式：（A-λI）X=0 （3.2）式（3.2）表示矩阵（A-λI）和非零向量X的乘积是零向量，式（3.2）有非零解的充分必要条件是矩阵（A-λI）是奇异的，即：det(A-λI)=0 该行列式可以表示为如下形式： a11–λa12 (1) a21 a22 –λ…a2n =0 (3.3) ………… A n1 a n2 …a nn 将式（3.3）中的行列式展开后，可以得到一个n阶多项式，称为特征多项式： f(λ)=det(A-λI)=(-1)n(λn+c1λn-1+c2λn-2+…+c n-1λ+c n) (3.4) n阶多项式一共有n个根（可以有重根），将每个根λ带入式（3.2），可以得到一个非零解向量。 习题：求下列矩阵的特征多项式的系数和特征值λj： 3 -1 0 A= -1 2 -1 0-1 3 解：在MATLAB中输入命令： A=【3 -1 0；-1 2 -1；0 -1 3】； c=poly(A) roots(c) 得到： 实验四：Lagrange插值多项式 实验目的：理解Lagrange插值多项式的基本概念，熟悉Lagrange插值多项式的公式源代码，并能根据所给条件求出Lagrange插值多项式，理解龙格现象。 %功能：对一组数据做Lagrange插值 %调用格式：yi=Lagran_(x,y,xi) %x,y：数组形式的数据表 %xi：待计算y值的横坐标数组 %yi：用Lagrange还擦之算出y值数组 function fi=Lagran_（x,f,xi） fi=zeros(size(xi)); np1=length(f); for i=1:np1 z=ones(size(xi)); for j=i:np1 if i~=j,z=z.*(xi-x(j))/(x(i)-x(j));end end fi=fi+z*f(i); end return 习题：已知4对数据（1.6，3.3），（2.7，1.22），（3.9，5.61），（5.6，2.94）。写出这四个数据点的Lagrange插值公式，并 1、 在MATLAB 中用Jacobi 迭代法讨论线性方程组， 1231231234748212515 x x x x x x x x x -+=?? -+=-??-++=? （1）给出Jacobi 迭代法的迭代方程，并判定Jacobi 迭代法求解此方程组是否收敛。 （2）若收敛，编程求解该线性方程组。 解（1）：A=[4 -1 1;4 -8 1;-2 1 5] %线性方程组系数矩阵 A = 4 -1 1 4 -8 1 -2 1 5 >> D=diag(diag(A)) D = 4 0 0 0 -8 0 0 0 5 >> L=-tril(A,-1) % A 的下三角矩阵 L = 0 0 0 -4 0 0 2 -1 0 >> U=-triu(A,1) % A 的上三角矩阵 U = 0 1 -1 0 0 -1 0 0 0 B=inv(D)*(L+U) % B 为雅可比迭代矩阵 B = 0 0.2500 -0.2500 0.5000 0 0.1250 0.4000 -0.2000 0 >> r=eigs(B,1) %B 的谱半径 r = 0.3347 < 1 Jacobi迭代法收敛。 (2)在matlab上编写程序如下： A=[4 -1 1;4 -8 1;-2 1 5]; >> b=[7 -21 15]'; >> x0=[0 0 0]'; >> [x,k]=jacobi(A,b,x0,1e-7) x = 2.0000 4.0000 3.0000 k = 17 附jacobi迭代法的matlab程序如下： function [x,k]=jacobi(A,b,x0,eps) % 采用Jacobi迭代法求Ax=b的解 % A为系数矩阵 % b为常数向量 % x0为迭代初始向量 % eps为解的精度控制 max1= 300; %默认最多迭代300，超过300次给出警告D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 B=D\(L+U); f=D\b; x=B*x0+f; k=1; %迭代次数 while norm(x-x0)>=eps x0=x; x=B*x0+f; k=k+1; if(k>=max1) disp('迭代超过300次，方程组可能不收敛'); return; end end MATLAB与数值分析实验报告 报告人：秦旸照 学号： 2015020901033 时间： 2016.4.8 电子科技大学电子工程学院 一、实验目的 实验一：MATLAB软件平台与程序设计实验 二、实验原理 1.熟练掌握矩阵的生成、加、减、乘、除、转置、行列式、逆、范数等运算操作。（用.m文件和Matlab函数编写一个对给定矩阵进行运算操作的程序） 2. 熟练掌握算术符号操作和基本运算操作，包括矩阵合并、向量合并、符号转换、展开符号表达式、符号因式分解、符号表达式的化简、代数方程的符号解析解、特征多项式、函数的反函数、函数计算器、微积分、常微分方程的符号解、符号函数的画图等。（用.m文件编写进行符号因式分解和函数求反的程序） 3. 掌握Matlab函数的编写规范。 4.掌握Matlab常用的绘图处理操作，包括：基本平面图、图形注释命令、三维曲线和面的填充、三维等高线等。（用.m文件编写在一个图形窗口上绘制正弦和余弦函数的图形，并给出充分的图形注释） 5. 熟练操作MATLAB软件平台，能利用M文件完成MATLAB的程序设计。 三、实验方案 1. 编程实现以下数列的图像，用户能输入不同的初始值以及系数。并以 x，y为坐标显示图像 x(n+1) = a*x(n)-b*(y(n)-x(n)^2); y(n+1) = b*x(n)+a*(y(n)-x(n)^2) 2. 编程实现奥运5环图，允许用户输入环的直径。 3. 实现对输入任意长度向量元素的冒泡排序的升序排列。 不允许使用sort函数。 四、实验结果 1. 编程实现以下数列的图像，用户能输入不同的初始值以及系数。并以 x，y为坐标显示图像 数值分析实验 ——MATLAB实现 姓名sumnat 学号2013326600000 班级13级应用数学2班 指导老师 2016年1月 一、插值：拉格朗日插值 (1) 1、代码： (1) 2、示例： (1) 二、函数逼近：最佳平方逼近 (2) 1、代码： (2) 2、示例： (2) 三、数值积分：非反常积分的Romberg算法 (3) 1、代码： (3) 2、示例： (4) 四、数值微分：5点法 (5) 1、代码： (5) 2、示例： (6) 五、常微分方程：四阶龙格库塔及Adams加速法 (6) 1、代码：四阶龙格库塔 (6) 2、示例： (7) 3、代码：Adams加速法 (7) 4、示例： (8) 六、方程求根：Aitken 迭代 (8) 1、代码： (8) 2、示例： (9) 七、线性方程组直接法：三角分解 (9) 1、代码： (9) 2、示例： (10) 八、线性方程组迭代法：Jacobi法及G-S法 (11) 1、代码：Jacobi法 (11) 2、示例： (12) 3、代码：G-S法 (12) 4、示例： (12) 九、矩阵的特征值及特征向量：幂法 (13) 1、代码： (13) 2、示例： (13) 一、插值：拉格朗日插值 1、代码： function z=LGIP(x,y)%拉格朗日插值 n=size(x); n=n(2);%计算点的个数 syms a; u=0;%拉格朗日多项式 f=1;%插值基函数 for i=1:n for j=1:n if j==i f=f; else f=f*(a-x(j))/(x(i)-x(j)); end end u=u+y(i)*f;f=1; end z=expand(u);%展开 2、示例： >> x=1:6; y1=x.^5+3*x.^2-6; y2=sin(x)+sqrt(x); >> f1=LGIP(x,y1) f1 = -6+3*a^2+a^5 %可知多项式吻合得很好 >> f2=vpa(LGIP(x,y2),3) f2 = .962e-1*a^4+1.38*a+.300*a^2+.504-.436*a^3-.616e-2*a^5数值分析MATLAB上机实验
matlab数值计算(命令与示例)
《MATLAB与数值分析》第一次上机实验报告
Matlab作业3(数值分析)答案
数值分析的matlab实现
matlab2012实验1参考答案
数值分析的MATLAB程序
同济大学数值分析matlab编程题汇编
数值分析实验报告1
实验6 Matlab数值计算实验报告
数值分析matlab代码
实验6答案 Matlab数值计算
第2讲 matlab的数值分析
MATLAB数值分析实验三(线性方程求解及精度分析)
数值分析幂法与反幂法-matlab程序
基于MATLAB数值分析实验报告
matlab数值分析例题
MATLAB与数值分析实验报告一
数值分析实验— MATLAB实现