数值分析MATLAB上机实验.docx

数值分析实习报告姓名：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 实验结果51.3.1 最正确平方逼近法51.3.2 拉格朗日插值法71.3.3 比照82 第二题92.1 实验目的92.2 实验原理和方法102.3 实验结果102.3.1 第一问102.3.2 第二问112.3.3 第三问113 第三题123.1 实验目的123.2 实验原理和方法123.3 实验结果124 MATLAB程序141第一题某过程涉及两变量x与y，拟分别用插值多项式与多项式拟合给出其对应规律的近似多项式，xi及yi之间的对应数据如下：12345678910x12345678910y请用次数分别为3，4，5，6的多项式拟合并给出最好近似结果f(x)。请用插值多项式给出最好近似结果。实验目的：学习逼近与插值的原理与编程方法，由给出的点构造多项式，在某个范围内近似代替点所代表的函数，以便于简化对未知函数的各种计算。试验原理与方法：实验原理：拉格朗日插值法中先构造插值根底函数：lkx=j=0jknx-xixk-xi k=0,1,2,n，然后构造出拉格朗日多项式：pnx=k=0nj=0jknx-xixk-xifxk。最正确平方逼近中，设逼近函数Pnx=a0+a1x+a2x2+anxn，逼近函数与真实函数之差r=Pnx-y，r1r2rn=11x1x2x1nx2n1xnxnna0a1an-y1y2yn，即：r=Xa-Y，根据最小二乘准那么令i=0nri2=min，可以得到a=XTX-1XTY。实验方法：逼近法采用最正确平方逼近，依据最小二乘原那么：i=0nri2=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.033*x3 + 19.33*x2 - 94.48*x拟合结果：12345678910x12345678910y四次多项式：- 0.3818*x4 + 7.368*x3 - 42.14*x2 + 73.53*x拟合结果：12345678910x12345678910y五次多项式：拟合结果：12345678910x12345678910y六次多项式：拟合结果：12345678910x12345678910y比照可知，六次多项式拟合结果最好。拉格朗日插值法插值多项式注：此多项式为拉格朗日多项式的近似式，当x=10的时候偏差可以到达23以上。比照数据：123456789xy1011121314151617xy插值结果：123456789xy1011121314151617xy其中红点表示参考点。比拟选取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,4T；b2=100,-200,345T。(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,-10T。(3)A行分别为A1=1,3,A2=-7,1；b1=4,6T。2.1试验目的学习jacobi迭代法与GuassSeidel迭代法的原理与编程方法，研究方程组系数阵与右边项对方程的解及其收敛性的影响，判断迭代法的收敛条件。实验原理与方法实验原理：将方程组系数阵A分解为A=D+L+U,其中D为对角阵，L为减去D的下三角阵，U为减去D的上三角阵。Jacobi迭代法中构造如下迭代公式：xk+1=-D-1L+Uxk+D-1b而Gauss-Seidel迭代法的迭代公式为：xk+1=-D+L-1Uxk+D+L-1b初始值直接选取为0。在判断其收敛性时，分别求解其迭代矩阵的谱半径G，G=max1inli， li为迭代矩阵的特征值。实验方法：分别编写jacobi迭代及其收敛判别函数与Seidel迭代及其收敛判别函数。如果在初试迭代步数之内还未收敛就进展收敛判别，收敛判别的依据是迭代矩阵的谱半径是否小于1。比拟同一方程组的jacobi迭代法与Seidel迭代法的结果是否一样，在到达精度要求后比拟两种方法的迭代次数，比拟哪一个的效率更高。比拟方程组系数阵与等号右边的变化会对方程的解与收敛速度造成什么影响。如果迭代不收敛，那么考虑为什么不收敛，如果把方程组系数阵进展强对角占优处理，是否会收敛。2.3实验结果规定误差界：1e-42.3.1第一问 A=62-11-341-24,b=-324.由jacobi迭代法求得x= -0.72730.80810.2525，设定迭代20次，实际迭代16次，精度为。由seidel迭代法求得x= -0.72720.80810.2525，设定迭代20次，实际迭代10次，精度为：A=62-11-341-24，b=100-200345由jacobi迭代法求得x= 36.3636-2.0707114.0404，设定迭代40次，实际迭代23次，精度为。由Seidel迭代法求得x= 36.3637-2.0707114.0404，设定迭代20次，实际迭代15次，精度为通过比照可知：1、Seidel迭代的收敛速度明显高于jacobi迭代。2、b矩阵对收敛速度与误差精度有影响，b中元素较大时会放慢收敛速度并加大误差。2.3.2第二问 A=10.80.80.80.810.80.81，b=321由jacobi迭代法求解，100次迭代尚且不能到达精度。此时调用jacobi迭代法的收敛判别函数，求得特征值为：l1=-1.6，l2、3=0.8，GJ=1.6>1，迭代不收敛。由于Seidel迭代法求解x= 5.76910.7693 -4.2307，迭代次数31，精度为 A=10.80.80.80.810.80.81，b=50-10Jacobi迭代不收敛。Seldel迭代法求得x= 32.69227.6922 -42.3076 ，迭代次数38，精度为。比拟得知A矩阵元素如果相差很小，迭代次数会大幅增加，综合比拟可知b矩阵元素如果相差很大会增加迭代次数。2.3.3第三问试验结果与讨论A=13-71，b=46此时GJ=4.5826>1，GG=21>1，jacobi迭代法与Seidel迭代法都不收敛。如果交换A中行的顺序，得到-7113，用jacobi迭代计算，迭代8次，解得x=-0.63641.5454。用Seidel迭代法计算，只需迭代5次，得到x= -0.6364 1.5455，精度为。此时GJ=0.2182>GG=0.0476，从此可以看出收敛速度的快慢。3第三题给定函数，及节点，求其三次样条插值多项式可取I型或II型边界条件，并画图及及的图形进展比拟分析。注：涉及到线性方程组求解问题需采用适当的求解算法。学习三弯矩法的原理与编程方法，比照原函数与三次样条插值的结果。实验原理与方法实验原理：记hi=xi-xi-1，i=hihi+hi+1，li=1-i，gi=6hi+hi+1yi+1-yihi+1-yi-yi-1hi，给出插值两端处的二阶导数S''x0=y''x0=M0 S''xn=y''xn=Mn。组成递推式：iMi-1+2Mi+liMi+1=gi由于系数阵按行对角占优，方程组存在唯一确定解，可以使用高斯列主元消去法来解方程。最后将各个参数带入样条函数Six=Mi-1xi-x36hi+Mix-xi-136hi+yi-1-Mi-16hi2xi-xhi+yi-Mi6hi2x-xi-1hi即可求得样条函数。实验方法由于在此题中x(i+1)-x(i)=1,所以h(i)=1。在编程中直接将h设置成常数，简化了运算。首先求解 m、l、g，然后列出方阵求解M(i)，在求解方程组的过程中采用列主元素高斯消去法，分为消元与回代两个过程，编写这两个函数，解出除了两端的M，而两端点的M值等于两端点的函数二阶导数值。编写函数求出样条函数的系数，然后求出方程，对于三弯矩法三次样条函数，如果有n个点，那么有n-1个样条函数，除了两端需要求解n-2个M值，即解n-2阶方程组。在表达样条函数的时候采用if语句，对不同的区间进展划分，然后细分-5,5这个区间，间隔将其分为100份，这样可以表达出连续性，此时绘图比照三次样条函数与原函数。此题中原函数的二阶导不是计算机解出的没有编写相关程序。此题中求解的样条函数，MATLAB系统自动的将公因式提取，并且合并同类项，所以表达出的函数并不整齐规律，为了更好地表达三次样条函数的构造与性质，我专门手写了规整的样条函数。实验结果三弯矩法求解代入x=-5 -4 -3 -2 -1 0 1 2 3 4 5y=fx=11+x2f''-5=f''50.00842y= 0 M= 求得样条函数系数阵为：S=样条函数：(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*X3)/5000 + 13119/10000(3199*(X - 1)3)/10000 - (9357*X)/10000 + (619*X3)/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)= *(-4-x)3+*(x+5)3+*(-4-x)+*(x+5) x-5,-4S2(x)= *(-3-x)3+*(x+4)3+*(-3-x)+*(x+4) x-4,-3S3(x)= *(-2-x)3+*(x+3)3+*(-2-x)+*(x+3) x-3,-2S4(x)= *(-1-x)3+*(x+2)3+*(-1-x)+*(x+2) x-2,-1S5(x)= *(-x)3-*(x+1)3+*(-x)+(x+1) x-1,0S6(x)= *(1-x)3+*x3+*(1-x)+*x x0,1S7(x)= *(2-x)3+*(x-1)3+*(2-x)+*(x-1) x1,2S8(x)= *(3-x)3+*(x-2)3+*(3-x)+*(x-2) x2,3S9(x)= *(4-x)3+*(x-3)3+*(4-x)+*(x-3) x3,4S10(x)= *(5-x)3+*(x-4)3+*(5-x)+*(x-4) x4,5输入：a=-5:0.1:5;for i=1:length(a)b(i)=f(a(i);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; endendb=0;for i=0:n for k=1:N b=b+y(k)*x(k)i end Q(i+1,1)=b; b=0;ends=inv(P)*Q;f=0;for i=1:n+1 f=f+s(i)*t(i-1); simplify(f);endf=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);endif(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<delta) break; end Pendx=x'k,errjacobi迭代收敛判别函数：function s=liansanxing(A)N,N=size(A);D=zeros(N);LU=zeros(N);for i=1:N D(i,i)=A(i,i);endLU=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<delta) break; endendPx=x'k,errseidel迭代收敛判别函数：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); endendfor i=1:N D(i,i)=A(i,i);endU=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);endfor 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);endA=zeros(n-2);for i=1:n-2 A(i,i)=2;endfor i=1:n-3 A(i,i+1)=v(i);endfor i=2:n-2 A(i,i-1)=u(i);endB=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);ends=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); endendx=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;endb(n-1)=mn/6;for i=1:n-2 b(i)=m(i)/6;endc(1)=y(1)-m0/6;for i=2:n-1 c(i)=y(i)-m(i-1)/6;endd(n-1)=y(n)-mn/6;for i=1:n-2 d(i)=y(i+1)-m(i)/6;endS=a b c d;S;for i=1:n-1s(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*x3+1.3119*(1-x)+0.3762*xelseif 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第 27 页