数值分析实验作业,gauss消去法的数值稳定性分析资料
Gauss列主元消去法实验
<数值计算方法>实验报告1.实验名称实验2 Gauss 列主元消去法2.实验题目用Gauss 列主元消去法求解线性方程组。
0.0011 2.0002 3.0003 1.0001.0001 3.7122 4.6233 2.0002.0001 1.0722 5.6433 3.000x x x x x x x x x ++=⎧⎪-++=⎨⎪-++=⎩3.实验目的加深自己对Gauss 列主元消去法的理解和认识，并且通过做实验或做练习来加强自己Gauss 列主元消去法的掌握，学会并灵活运用Gauss 列主元消去法来求解方程组。
4.基础理论-------Gauss 列主元消去法1.Gauss 列主元消去法的基本思想是：在进行第k （k=1,2，...，n-1)步消元时，从第k 列的kk a 及以下的各元素中选取绝对值最大的元素，然后通过行变换将它交换到主元素kk a 的位置上，再进行消元。
2.Gauss 列主元消去法的优点：当kk a （k=1,2，...，n-1)的绝对值很小时，用Gauss 列主元消去法来求解方程组时，可以避免所的数值结果产生较大误差或失真。
5.实验环境实验系统：Win 7实验平台：VisualC++语言6.实验过程写出算法→编写程序→计算结果Gauss 列元消去法的算法Input:方程组未知量的个数n;增广矩阵()()1,2,...,T ij A a A A An ==，其中i=1,2,…,n; j=1,2,…,n+1Output:方程组的解x1,x2,…,xn,或失败信息。
1. for i ←1ton-1 do;2. temp ←|ii a |;3. p ←I;4. for j ←i+1 to n do5. if ||ji a >temp then6. p ←j;8. end9. end10. if temp=0 then11. |return False;12. end13. if p ≠I then14. p A ⇔i A ;//i,p 两行交换15. end//列选主元16. for j ←i+1 to n do17.*j ji i A m A -ji m ←/ji ii a a ;18. j A ←*j ji i A m A -;//消元19. end7.实验结果原方程组的解为：X1=-0.490396 , x2=-0.051035 ,x3=0.3675208.附录程序清单#include<iostream.h> #include"stdio.h"#include"math.h"void main ( ){ int n=3,i,j,k,p;doubleA[10][10]={{0.001,2.000,3.000,1.000},{-1.000,3.712,4.623,2.000},{-2.0 00,1.072,5.643,3.000}},temp,m,x[100];for(i=0;i<n;i++){ //选主元temp=fabs(A[i][i]); p=i;for(k=i+1;k<n;k++)if(fabs(A[k][i])>temp){temp=fabs(A[k][i]); p=k;}if(temp==0){ printf("\n无法求解：");return;}if(p!=i)for(j=0;j<n+1;j++){ temp=A[i][j];A[i][j]=A[p][j];A[p][j]=temp;}//消元for(k=i+1;k<n;k++){ m=A[k][i]/A[i][i];for(j=i+1;j<=n;j++)A[k][j]=A[k][j]-m*A[i][j];}}//回代for(i=n-1;i>=0;i--){x[i]=A[i][n];for(j=i+1;j<n;j++)x[i]=x[i]-A[i][j]*x[j];x[i]=x[i]/A[i][i];}printf("\nx=\n");for(i=0;i<n;i++)printf("%f \n",x[i]);}。
数值分析实验报告
天津工业大学数值分析实验报告电子与信息工程学院例题一用Gauss 消去法来解方程组x 1+x 2+x 3=612x 1−3x 2+3x 3=15−18x 1+3x 2−x 3=−15模型：高斯（Gauss ）消去法高斯消去法的基本思想是：先逐次消去变量，将方程组化成同解的上三角形方程组，此过程成为消元过程。
然后按方程相反顺序求解上三角方程组，得到原方程组的解，此过程称为回代过程。
这种方法称为Gauss 消去法。
将方程组改成如下形式a 11 1x 1+a 12 1x 2+······+a 1n 1xn =b 11a 21 1 x 1+a 22 1 x 2+······+a 2n 1 xn =b 2 1 : ∶ ∶ ∶ : ∶ ∶ ∶ a n 1 1x 1+a n 2 1 x 2+······+a nn 1xn =b n1 （2-1） 消元过程：第一步：设a 11 1≠0，记l i 1=a i 1/a 11 1(i =2,3,······，n )，方程中第i 个方程减去第一个方程乘以l i 1，完成第一次消元。
将上述方程式可以改写为a 111x 1+a 12 1x 2+······+a 1n 1xn =b 1 1a 22 2 x 2+······+a 2n 2 xn =b 2 2 ∶ ∶ ∶ ∶ ∶ ∶ a n 22 x 2+······+a nn 2 xn =b n 2（2-2） 其中a ij (2)=a ij (1)−l i 1a 1j (1)，b i (2)=b i(1)−l i 1b 1 1(i ,j =2,3·········，n )。
数值分析实验作业-gauss消去法的数值稳定性分析说课讲解
数值分析实验作业-g a u s s消去法的数值稳定性分析实验3.1 Gauss 消去法的数值稳定性试验实验目的:观察和理解Gauss 消元过程中出现小主元(即)(k kka 很小)时引起的方程组解的数值不稳定性。
实验内容：求解方程组b Ax =，其中（1）⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⨯=11212592.1121-130.6-291.51314.59103.015-1A ，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=2178.4617.591b ； （2）⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----=2010151526990999999999.23107102A ,⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=15019000000000.582b .实验要求：（1） 计算矩阵的条件数，判断系数矩阵是良态的还是病态的。
（2） 用Gauss 列主元消去法求得L 和U 及解向量421,R x x ∈.（3） 用不选主元的Gauss 消去法求得L ~和U ~及解向量421~,~R x x ∈.（4） 观察小主元并分析其对计算结果的影响. 程序如下:计算矩阵条件数及Gauss 列主元消去法: format longengA1=[0.3e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1];b1=[59.17;46.78;1;2]; n=4;k2=cond(A1) %k2为矩阵的条件数；for k=1:n-1a=max(abs(A1(k:n,k))); [p,k]=find(A1==a); B=A1(k,:);c=b1(k);A1(k,:)=A1(p,:);b1(k)=b1(p); A1(p,:)=B;b1(p)=c; if A1(k,k)~=0A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k);A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n); else break end endL1=tril(A1,0); for i=1:n L1(i,i)=1; end L=L1U=triu(A1,0) for j=1:n-1b1(j)=b1(j)/L(j,j);b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j); endb1(n)=b1(n)/L(n,n); for j=n:-1:2b1(j)=b1(j)/U(j,j);b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j); endb1(1)=b1(1)/U(1,1); x1=b1运行结果如下： K2=68.43；⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡---⨯=-14929.00202.00893.0011755.04724.00011079.2600118L⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=801.0000231.1835.2001314.5902592.11U 1x =[18.9882;3.3378;-34.747;-33.9865] 不选主元的Gauss 消去法程序：clearformat longengA1=[0.3e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1]; b1=[59.17;46.78;1;2]; n=4;for k=1:n-1A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k);A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n); endL1=tril(A1,0); for i=1:n L1(i,i)=1; end L=L1U=triu(A1,0) for j=1:n-1b1(j)=b1(j)/L(j,j);b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j); endb1(n)=b1(n)/L(n,n); for j=n:-1:2b1(j)=b1(j)/U(j,j);b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j); endb1(1)=b1(1)/U(1,1); x1=b1程序运行结果如下：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⨯⨯=10189.010333.3011168.21033.3700110637.170001~151515L⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡--⨯-⨯-⨯-⨯=-5.000816010637.171091.5210043.101314.59103.0~15151815U ]0;0;0005.1;6848.23[~1=x同理可得2A 对应的系数矩阵条件数及Gauss 列主元消去法求解结果： K2=8.994;⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-⨯=1333.04.0001104.0-3.0-0015.0000112-L ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=36667.300030.26005.155.2010710U ]0.1;;0.1;0.1;10444.0[152-⨯=-x不选主元的Gauss 消去法结果：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⨯=1400.0109999.0-001104998.2-5.00013.0-0001~1212L ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⨯⨯--=-3667.3000107495.5109987.14003.26100.1010710~121212U ]000145.1;99994.0;000.1;1045.1[~52-⨯-=-x实验4.5 三次样条插值函数的收敛性问题提出：多项式插值不一定收敛的，即插值的节点多，效果不一定就好。
高斯消去法的数值稳定性实验报告
end
%½«³£ÊýÏî´ÓAÖзÖÀë
X=A(:,length(A(1,:)));
A=A(:,1:length(A(1,:))-1);
%µÃµ½LÓëU¾ØÕó
fori=1:length(A(:,1))
forj=1:length(A(i,:))
if(i==j)
L(i,j)=1;
U(i,j)=A(i,j);
00-2.83542034035461-0.439277018213440
0000.724838760535495
X2=
5.90739256961544
1.88284104193886
-22.2540722833854
14.5809976298923
(3)beselect.m文件与noselect.m文件
%½«³£ÊýÏî¼Ó½øÈ¥
A(:,length(A(1,:))+1)=Y;
%´ÓÉÏÍùÏÂÑ »·
fori=1:length(A(:,1))
%Ô¤¼ÆËãÿÐеÚÒ»¸öÊý×Ö
fork=i:length(A(:,i))
form=1:i-1
A(k,i) = A(k,i)-A(k,m)*A(m,i);
end
end
end
%¼ÆËãXÖµ
fori=length(X):-1:1
forj=i+1:length(U(i,:))
X(i)=X(i)-U(i,j)*X(j);
end
X(i)=X(i)/U(i,i);
end
end
noselect.m文件：
function[ L,U,X ] = noselect( A,Y )
数值分析实验6 (1)
Lab06．Gauss 列主元素消去法实验【实验目的和要求】1．使学生深入理解并掌握Gauss 消去法和Gauss 列主元素消去法步骤；2．通过对Gauss 消去法和Gauss 列主元素消去法的程序设计，以提高学生程序设计的能力；3．对具体问题，分别用Gauss 消去法和Gauss 列主元素消去法求解。
通过对结果的分析比较，使学生感受Gauss 列主元素消去法优点。
【实验内容】1．根据Matlab 语言特点，描述Gauss 消去法和Gauss 列主元素消去法步骤。
2．编写用不选主元的直接三角分解法解线性方程组Ax=b 的M 文件。
要求输出Ax=b 中矩阵A 及向量b ，A=LU 分解的L 与U ，det A 及解向量x 。
3．编写用Gauss 列主元素消去法解线性方程组Ax=b 的M 文件。
要求输出Ax=b 中矩阵A 及向量b 、PA=LU 分解的L 与U 、det A 及解向量x ，交换顺序。
4．给定方程组(1) ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--11134.981.4987.023.116.427.199.103.601.3321x x x (2) ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----15900001.582012151********.23107104321x x x x 先用编写的程序计算，再将(1)中的系数3.01改为3.00，0.987改为0.990；将(2)中的系数2.099999改为2.1，5.900001改为9.5,再用Gauss 列主元素消去法解，并将两次计算的结果进行比较。
【实验仪器与软件】1．CPU 主频在1GHz 以上，内存在128Mb 以上的PC ；2．Matlab 6.0及以上版本。
实验讲评：实验成绩：评阅教师：2011年6 月 30日Lab06．Gauss列主元素消去法实验一算法描述1编写用不选主元的直接三角分解法解线性方程组Ax=b的M文件程序如下function [x,l,u]=malu(A,b)format shortn=length(b);u=zeros(n,n);l=eye(n,n);u(1,:)=A(1,:);l(2:n,1)=A(2:n,1)/u(1,1);for k=2:nu(k,k:n)=A(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n);l(k+1:n,k)=(A(k+1:n,k)-l(k+1:n,1:k-1)*u(1:k-1,k))/u(k,k);endy=zeros(n,1);y(1)=b(1);for k=2:ny(k)=b(k)-l(k,1:k-1)*y(1:k-1);endx=zeros(n,1);x(n)=y(n)/u(n,n);for k=n-1:-1:1x(k)=(y(k)-u(k,k+1:n)*x(k+1:n))/u(k,k);end2、Gauss列主元消去法解线性方程组程序如下：function [Determ,x]=magauss2(A,b,flag)%Gauss列主元素消去法解线性方程组Ax=b，A为系数矩阵，b为右端项%若flag=0，不显示中间消去过程，否则显示中间消去过程，默认为0 %输出项Determ为矩阵A的行列式值，x为解向量if nargin<3,flag=0;endDeterm=1;n=length(b);for k=1:(n-1)[ap,p]=max(abs(A(k:n,k)));p=p+k-1;if ap==0printf('divide by zero!');Determ=0;end%换行if p>kt=A(k,:); A(k,:)=A(p,:);A(p,:)=t;t=b(k);b(k)=b(p);b(p)=t;Determ=-Determ;end%消元计算m=A(k+1:n,k)./A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n); b(k+1:n)=b(k+1:n)-m*b(k);A(k+1:n,k)=zeros(n-k,1);if flag~=0,Ab=[A,b],end %展示消元过程Determ=A(k,k).* Determ;endif A(n,n)==0printf('divide by zero!');Determ=0;end%回代求解x=zeros(n,1);x(n)=b(n)/A(n,n);for i=(n-1):-1:1x(i)=(b(i)-A(i,i+1:n)*x(i+1:n))/A(i,i);endDeterm=A(n,n).* Determ;三计算过程;直接三角分解法(1):A=[3.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34]b=[1;1;1][x1,l1,u1]=malu(A,b);x1l1u1A =3.0100 6.0300 1.99001.2700 4.1600 -1.23000.9870 -4.8100 9.3400b =111x1 =1.0e+003 *1.5926-0.6319-0.4936l1 =1.0000 0 00.4219 1.0000 00.3279 -4.2006 1.0000u1 =3.0100 6.0300 1.99000 1.6158 -2.06960 0 -0.0063>>A=[3.00 6.03 1.99;1.27 4.16 -1.23;0.990 -4.81 9.34] b=[1;1;1][x1,l1,u1]=malu(A,b);x1l1u1A =3.0000 6.0300 1.99001.2700 4.1600 -1.23000.9900 -4.8100 9.3400b =111x1 =119.5273-47.1426-36.8403l1 =1.0000 0 00.4233 1.0000 00.3300 -4.2306 1.0000u1 =3.0000 6.0300 1.99000 1.6073 -2.07240 0 -0.0844>> (2)A=[10,-7,0,1;-3,2.09999,6,2;5,-1,5,-1;2,1,0,2] b=[8;5.900001;5;1][x1,l1,u1]=malu(A,b);x1l1u1A =10.0000 -7.0000 0 1.0000-3.0000 2.1000 6.0000 2.00005.0000 -1.0000 5.0000 -1.00002.0000 1.0000 0 2.0000b =8.00005.90005.00001.0000x1 =0.0000-1.00001.00001.0000l1 =1.0e+005 *0.0000 0 0 0-0.0000 0.0000 0 00.0000 -2.5000 0.0000 00.0000 -2.4000 0.0000 0.0000u1 =1.0e+006 *0.0000 -0.0000 0 0.00000 -0.0000 0.0000 0.00000 0 1.5000 0.57500 0 0 0.0000>> A=[10,-7,0,1;-3,2.1,6,2;5,-1,5,-1;2,1,0,2]b=[8;5.9;5;1][x1,l1,u1]=malu(A,b);x1l1u1A =10.0000 -7.0000 0 1.0000-3.0000 2.1000 6.0000 2.00005.0000 -1.0000 5.0000 -1.00002.0000 1.0000 0 2.0000b =8.00005.90005.00001.0000x1 =NaNNaNNaNNaNl1 =1.0000 0 0 0-0.3000 1.0000 0 00.5000 Inf 1.0000 00.2000 Inf NaN 1.0000u1 =10.0000 -7.0000 0 1.00000 0 6.0000 2.30000 0 -Inf -Inf0 0 0 NaN高斯列主元消去法：（1）：A1=[3.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34] b=[1;1;1]condA1=cond(A1)%(1)中A的条件数[Determ1,x1]=magauss2(A1,b,1)A2=[3.00 6.03 1.99;1.27 4.16 -1.23;0.990 -4.81 9.34] b=[1;1;1]condA2=cond(A2)[Determ2,x2]=magauss2(A2,b,1 )运行结果如下A1 =3.0100 6.0300 1.99001.2700 4.1600 -1.23000.9870 -4.8100 9.3400b =111condA1 =3.0697e+004Ab =3.0100 6.0300 1.9900 1.00000 1.6158 -2.0696 0.57810 -6.7873 8.6875 0.6721 Ab =3.0100 6.0300 1.9900 1.00000 -6.7873 8.6875 0.67210 0 -0.0015 0.7381 Determ1 =-0.0305x1 =1.0e+003 *1.5926-0.6319-0.4936A2 =3.0000 6.0300 1.99001.2700 4.1600 -1.23000.9900 -4.8100 9.3400b =111condA2 =2.3028e+003Ab =3.0000 6.0300 1.9900 1.00000 1.6073 -2.0724 0.57670 -6.7999 8.6833 0.6700Ab =3.0000 6.0300 1.9900 1.00000 -6.7999 8.6833 0.67000 0 -0.0200 0.7350Determ2 =-0.4070x2 =119.5273-47.1426-36.8403（2）A1=[10,-7,0,1;-3,2.09999,6,2;5,-1,5,-1;2,1,0,2] b=[8;5.900001;5;1]condA1=cond(A1)%(1)中A的条件数[Determ1,x1]=magauss2(A1,b,1)A2=[10,-7,0,1;-3,2.1,6,2;5,-1,5,-1;2,1,0,2]b=[8;5.9;5;1]condA2=cond(A2)-[Determ2,x2]=magauss2(A2,b,1)A1 =10.0000 -7.0000 0 1.0000-3.0000 2.1000 6.0000 2.00005.0000 -1.0000 5.0000 -1.00002.0000 1.0000 0 2.0000b =8.00005.90005.00001.0000condA1 =5.8745Ab =10.0000 -7.0000 0 1.0000 8.00000 -0.0000 6.0000 2.3000 8.30000 2.5000 5.0000 -1.5000 1.00000 2.4000 0 1.8000 -0.6000 Ab =10.0000 -7.0000 0 1.0000 8.00000 2.5000 5.0000 -1.5000 1.00000 0 6.0000 2.3000 8.30000 0 -4.8000 3.2400 -1.5600 Ab =10.0000 -7.0000 0 1.0000 8.00000 2.5000 5.0000 -1.5000 1.00000 0 6.0000 2.3000 8.30000 0 0 5.0800 5.0800 Determ1 =-762.0009x1 =0.0000-1.00001.00001.0000A2 =10.0000 -7.0000 0 1.0000-3.0000 2.1000 6.0000 2.00005.0000 -1.0000 5.0000 -1.00002.0000 1.0000 0 2.0000b =8.00005.90005.00001.0000condA2 =5.8746Ab =10.0000 -7.0000 0 1.0000 8.00000 0 6.0000 2.3000 8.30000 2.5000 5.0000 -1.5000 1.00000 2.4000 0 1.8000 -0.6000 Ab =10.0000 -7.0000 0 1.0000 8.00000 2.5000 5.0000 -1.5000 1.00000 0 6.0000 2.3000 8.30000 0 -4.8000 3.2400 -1.5600 Ab =10.0000 -7.0000 0 1.0000 8.00000 2.5000 5.0000 -1.5000 1.00000 0 6.0000 2.3000 8.30000 0 0 5.0800 5.0800Determ2 =-762x2 =0.0000-1.00001.00001.0000四算法分析在（1）中；直接三角解法与高斯列主元消去解法的结果是一样的，而经过改变某些值后，其结果还是一样，看不出来那个方法好在（2）中直接三角解法出现结果：五总结通过实验，深入并基本掌握高斯消去法和高斯列主元消去法的步骤，二经过实际问题的求解，感受到高斯列主元消去法的优点，以及算法的不同产生的误差很大，对我们的问题求解有很大的影响。
计算方法实验报告(1)----Gauss消去法
计算方法实验报告实验名称：实验（一）Gauss 消去法班级：学生姓名：学号：班级序号：课内序号:指导老师：2018-2019学年第2学期一、实验名称：Gauss消去法二、实验学时: 2学时三、实验目的和要求1、掌握高斯消去法基础原理2、掌握高斯消去法解方程组的步骤3、能用程序语言对Gauss消去法进行编程实现四、实验过程代码及结果1、代码：using System;using System.Collections.Generic;using System.Linq;using System.Text;namespace ConsoleApplication_Gauss{class Program{//回带求值的过程static void CalcX(double[,] a, double[] x, int n){for (int i = n - 1; i >= 0; i--){double sum = 0;for (int j = i + 1; j < n; j++){sum += a[i, j] * x[j];}x[i] = (a[i, n] - sum) / a[i, i];}}//消元的过程static void CalcA(double[,] a, int n){for (int k = 0; k < n - 1; k++){for (int i = k + 1; i < n; i++){//double Lik = a[i, k] / a[k, k];// for (int j = k ; j <= n; j++)for (int j = n; j >= k; j--){a[i, j] = a[i, j] - a[i, k] / a[k, k] * a[k, j];}//a[i, k] = 0;}//Output}}//输出未知数x的值static void Output(double[] x, int n){for (int i = 0; i < n; i++){Console.WriteLine("x[{0}]={1}", i, x[i]);}}static void Output(double[,] a, int n){for (int i = 0; i < n; i++){//string s="";for (int j = 0; j <= n; j++){//s += string.Format("{0,-4}", a[i, j]);Console.Write("{0,6}", a[i, j]);}Console.WriteLine();}}//输入函数，表示输入一串值作为方程组的系数static void Input(double[,] a, int n){for (int i = 0; i <= n - 1; i++){string s = Console.ReadLine();string[] ss = s.Split(' ');for (int j = 0; j <= n; j++){a[i, j] = Convert.ToDouble(ss[j]);}}}static void Main(string[] args){Console.WriteLine("请输入矩阵的维数:");int n =Convert.ToInt32( Console.ReadLine());double[,] a = new double[n,n+1];Console.WriteLine("请输入矩阵的各个元;");Input(a, n);Console.WriteLine("------A（i,j)----------");Output(a, n);CalcA(a, n);Console.WriteLine("------消元之后A（i,j)----------");Output(a, n);double[] x = new double[n];CalcX(a, x, n);Output(x, n);Console.ReadLine();}}}2、结果：…。
数值计算方法课后实验-李维国
实验3.1 Gauss消去法的数值稳定性实验实验目的：观察和理解高斯消元过程中出现小主元(即|a(k)kk|很小)时引起方程组解数值不稳定性．实验内容：求解方程组Ax=b,其中实验要求：（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的．（2）用高斯列主元消去法求得L和U及解向量x1,x2 ．（3）用不选主元的高斯消去法求得L和U及解向量x1,x2 ．（4）观察小主元并分析对计算结果的影响．解：（1）cond(A1)=68.4296cond(A2)=8.9939A1矩阵条件数远大于1，A1矩阵为病态;A2矩阵条件数没有远远大于1，A2矩阵为良态。
（2）高斯列主元消去法：程序如下ClearA1=[0.3*power(10,-15),59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1];b1=[59.17;46.18;1;2];% A1=[1,2,-1,1;1,1,2,-1;3,-1,1,1;2,1,3,-1];%b1=[1;-2;6;-1];A2=[10,-7,0,1;-3,2.0999********,6,2;5,-1,5,-1;0,1,0,2];b2=[8;5.900000000001;5;1];A=A1; b=b1;n=input('n=');for k=1:n-1 [a,b3]=max(A(k:n,k)); A([k k-1+b3],:)=A([k-1+b3 k],:); b([k k-1+b3],:)=b([k-1+b3 k],:); if A(k,k)~=0A(k+1:n,k)=A(k+1:n,k)./A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); else stop end endA; b; L=tril(A,-1);U=triu(A,0); for i=1:1:n L(i,i)=1;end L; for j=1:n-1 b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);end b(n)=b(n)/L(n,n); b; y=b;for j=n:-1:2 y(j)=y(j)/U(j,j); y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);end y(1)=y(1)/U(1,1);结果如下:方程组一：L1： U1：解向量x1结果如下：方程组二结果如下：L2：U2:解向量x2结果如下：（3）不选主元的分解程序如下：ClearA1=[0.3*power(10,-15),59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1]; b1=[59.17;46.18;1;2];%A1=[6,2,1,-1;2,4,1,0;1,1,4,-1;-1,0,-1,3];%b1=[6;-1;5;-5];A2=[10,-7,0,1;-3,2.0999********,6,2;5,-1,5,-1;0,1,0,2];b2=[ 8;5.900000000001;5;1]; A=A2; b=b2; n=input('n='); for k=1:n-1A(k+1:n,k)= A(k+1:n,k)/A(k,k); A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); end L=tril(A,-1); U=triu(A,0); for i=1:1:n L(i,i)=1; end L ;for j=1:n-1 b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j); end b(n)=b(n)/L(n,n); b; y=b;for j=n:-1:2 y(j)=y(j)/U(j,j); y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j); end y(1)=y(1)/U(1,1);方程组一结果：L1: U1:解向量x1结果如下：方程组二结果：L2: U2:解向量x2结果如下：(4) 观察方程在两种不同方法下的结果可知：由于计算机字长是一定的，小主元会造成大数除以小数的结果超出字长，结果发生很大的变化。
数值分析实验报告---高斯消去法 LU分解法
数值分析实验报告---高斯消去法 LU分解法实验一：高斯消去法一、实验目的1. 掌握高斯消去法的原理2. 用高斯消去法解线性方程组3. 分析误差二、实验原理高斯消去法(又称为高斯-约旦消去法)是一种利用矩阵消元的方法，将线性方程组化为改进的阶梯形式，从而解出线性方程组的解的方法。
具体而言，高斯消去法将线性方程组的系数矩阵化为一个上三角矩阵，再利用回带法求解线性方程组的解。
三、实验内容1.1、用高斯消去法解线性方程组在具体实验中，我们将使用高斯消去法来解决下述的线性方程组。
5x+2y+z=102x+6y+2z=14x-y+10z=25为了使用高斯消去法来解这个方程组，首先需要将系数矩阵A进行变换，消除A矩阵中第一列中的下角元素，如下所示:1, 2/5, 1/50, 28/5, 18/50, 0, 49/28接着使用回代法来计算该方程组的解。
回代法的过程是从下往上进行的，具体步骤如下：第三个方程的解：z=49/28；第二个方程的解： y=（14-2z-2x）/6；第一个方程的解： x=（10-2y-z）/5。
1.2、分析误差在使用高斯消去法求解线性方程组时，一般会出现截断误差，导致得到的解与真实解之间存在一些误差。
截断误差的大小和矩阵的维数有关。
为了估计截断误差，我们使用矩阵B来生成误差，在具体实验中，我们将使用下面的矩阵：我们来计算该矩阵的行列式，如果方程组有唯一解，则行列性不为0。
本例中，行列式的值是 -1，因此方程组有唯一解。
然后我们计算真实解和高斯消去法得到的解之间的误差，具体公式如下所示：误差 = 真实解的范数 - 高斯消去法得到的解的范数其中，范数的定义如下：||x||1=max{|xi|}; ||x||2=sqrt{(|x1|^2 + |x2|^2 + ... + |xn|^2)}四、实验步骤1、将高斯消去法的每一个步骤翻译成代码，并保存为一个独立的函数。
2、将代码上传至 Python 交互式环境，并使用高斯消去法来解线性方程组。
高斯消去实验报告分析
一、实验背景高斯消去法（Gaussian Elimination）是线性代数中解决线性方程组的重要方法之一。
它通过初等行变换将系数矩阵化为行阶梯形矩阵，从而求解线性方程组。
本实验旨在通过编写程序实现顺序高斯消去法和列主元高斯消去法，并对比分析两种方法的计算效率和精度。
二、实验目的1. 理解高斯消去法的原理和步骤。
2. 编写程序实现顺序高斯消去法和列主元高斯消去法。
3. 比较两种方法的计算效率和精度。
4. 分析高斯消去法的适用范围和局限性。
三、实验内容1. 顺序高斯消去法- 实现顺序高斯消去法，将系数矩阵化为行阶梯形矩阵。
- 求解线性方程组，得到未知数的解。
2. 列主元高斯消去法- 实现列主元高斯消去法，将系数矩阵化为行阶梯形矩阵。
- 求解线性方程组，得到未知数的解。
3. 比较两种方法- 对比两种方法的计算效率和精度。
- 分析两种方法的优缺点。
四、实验步骤1. 编写程序- 使用C/C++、Python等编程语言编写顺序高斯消去法和列主元高斯消去法程序。
- 程序应包含以下功能：- 输入线性方程组的系数矩阵和常数项。
- 实现顺序高斯消去法和列主元高斯消去法。
- 输出求解结果。
2. 测试程序- 使用一组已知线性方程组测试程序的正确性。
- 验证求解结果与精确解是否一致。
3. 比较分析- 使用不同的线性方程组测试两种方法的计算效率和精度。
- 比较两种方法的运行时间、解的精度和误差。
五、实验结果与分析1. 顺序高斯消去法- 计算效率：顺序高斯消去法的时间复杂度为O(n^3)，其中n为方程组中未知数的个数。
- 精度：顺序高斯消去法可能会受到舍入误差的影响，导致求解结果与精确解存在一定的误差。
2. 列主元高斯消去法- 计算效率：列主元高斯消去法的时间复杂度同样为O(n^3)，但通过选择主元，可以减少舍入误差，提高求解精度。
- 精度：列主元高斯消去法具有更高的精度，求解结果与精确解更接近。
3. 比较分析- 顺序高斯消去法和列主元高斯消去法在计算效率上基本相同，但列主元高斯消去法具有更高的精度。
数值分析实验二（列主元Gauss消去法）
《数值分析》实验报告实验编号：实验二课题名称：列主元Gauss消去法一、算法介绍1、输入矩阵的阶数n，方程组的增广矩阵A；2、对k=0,1,…,n-2，循环：选取列中绝对值最大的元素，将主元所在的行的元素保存在数组temp[n+1]中。
若主元为零，则系数矩阵奇异，计算停止；否则，顺序进行。
如果绝对值最大的元素就在矩阵的对角线上，则进行普通高斯消元法的第一大步，否则将方程组系数换行之后再进行普通高斯消元法的第一大步；3、然后利用回代法求解线性方程组。
二、程序代码#include<iostream>#include<cmath>#include<iomanip>using namespace std;int main(){int n=0,k=0,i=0,j=0,h=0,g=0,flag=0,i1,j1;double max=0,m=0;cout<<"***利用列主元Gauss消元法求解线性方程组***"<<endl;cout<<"请输入矩阵的阶数："<<endl;cin>>n;double a[n][n+1];double t[n+1];double x[n];memset(a,0,sizeof(a));memset(x,0,sizeof(x));cout<<"请输入方程组的增广矩阵："<<endl;for(i=0;i<n;i++){for(j=0;j<n+1;j++){cin>>a[i][j];}}for(k=0;k<n-1;k++){max=0;j1=0;for(i=k;i<n;i++){if(fabs(a[i][k])>max){max=fabs(a[i][k]);i1=i;j1=k;}}if(max==0){cout<<"该系数矩阵为奇异矩阵，计算停止"<<endl;flag=1;break;}else{cout<<"第"<<j1+1<<"列中绝对值最大的元素是"<<a[i1][j1]<<"，在线性方程组的第"<<i1+1<<"行"<<endl;if(i1!=k){for(j=0;j<=n;j++){t[j]=a[i1][j];a[i1][j]=a[k][j];a[k][j]=t[j];}}for(i=k+1;i<=n-1;i++){m=a[i][k]/a[k][k];for(j=k;j<=n;j++)a[i][j]=a[i][j]-m*a[k][j];for(g=0;g<n;g++){for(h=0;h<n+1;h++)cout<<setiosflags(ios::fixed)<<setprecision(2)<<a[g][h]<<" ";cout<<endl;}cout<<endl;}}}if(flag==0){x[n-1]=a[n-1][n]/a[n-1][n-1] ;double sum=0;for(k=n-2;k>=0;k--){sum=0;for(i=n-1;i>=k;i--)sum+=a[k][i]*x[i];x[k]=(a[k][n]-sum)/a[k][k];}cout<<"该线性方程组的解为："<<endl;for(i=0;i<n;i++)cout<<"x"<<i+1<<"="<<setiosflags(ios::fixed)<<setprecision(2)<<x[i]<<endl;}system("pause");return 0;}三、运算结果截屏四、算法分析列主元Gauss消元法避免了普通高斯消元法中出现的问题：遇到某个主元为零或者当主元绝对值很小时，计算将会停止或求出的结果将与其实际结果相差很远。
数值分析A实验报告
的解。A 及 b 分别为由函数 matrix31 及 array31 所生成的方阵与向量。
function[x,det,flag]=gaussauto(A,b)
[n,m]=size(A);nb=length(b);
if n~=m error('The rows and columns of matrix A must be equal!');
结果分析，由于最大的主元与最小的主元的绝对值没有过于接近于 0，都 是“比较好的一个数”，然而，如果选择 0 能够作为模最小的数的话，则计算 过程将会崩溃。由此可见，虽然矩阵 A 的条件数非常大即显得相当“病态”， 然而，假设不使用 0 元素作为主元素的话，不管以最大模或者最小模作为主元
7
素的话，对于结果都是没有影响的，原因在于不管是哪里一种方法，选取的主 元的绝对值与其他元素相较起来都属于“良态”，当然这个是不考虑 0 元素的 存在的情况。
1.5.4 选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数，重复上述实 验，观察记录并分析实验的结果。
8
我们这里直接引用实验 2 即将出现的 Hilbert 矩阵作为系数矩阵，并给定 其有精确解[1,…,1]’。
我们考察六维度的 Hilbert 系数矩阵，有如下表达，即 H 为
设有
1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909
高斯消元法_实验报告
- - 华中科技大学数值分析实验报告系、年级研究生院2012级****类别硕士2013年5月6日实验6.1实验要求：根据教材实验6.1做出相应改编：分别使用Gauss 消元、列选主元。
全选主元的方法求解线性方程组，分别比拟三种消元方法的结果和算法的区别，并说明主元的选取在Gauss 消元的中的作用。
问题提出：Gauss 消去法是我们在线性代数中已经熟悉的。
但由于计算机的数值运算是在一个有限的浮点数集合上进展的，如何才能确保Gauss 消去法作为数值算法的稳定性呢？Gauss 消去法从理论算法到数值算法，其关键是主元的选择。
主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。
一般来说书本上采用的列选主元的方法对其线性方程组进展求解的，则我们是否可以选择一种行列都选取主元消去的方法来减小相应的误差呢？全主元消元法和列主元消元法一样都是由高斯消元法演变而来。
只不过选取主元的*围有所加大。
全选主元相对于列选主元的更加复杂化了，因为在运算的过程中导致了元的位置发生了变化，这样我们就不得不追踪每个元的位置。
本次实验就几个问题进展了matlab 实验分析，比拟几种计算方法的优劣性。
实验内容：考虑线性方程组编制一个程序：分别能进展Gauss 消去、列选主元Gauss 消去、全选主元Gauss 消去法进展解线性方程组。
对三种算法所得到的结果进展比拟，分析三种计算方法的准确性。
具体内容：〔1〕取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816 b A ，则方程有解T x )1,,1,1(* =。
取n=10、n=20计算矩阵的条件数。
分别编写利用matlab 编写运算程序，实现Gauss消去、列选主元消去以及全选主元消去的方法。
比拟三种计算方法的运算结果。
在列选主元的过程中分别采用每步消去过程总选取按模最小或按模尽可能小的元素作为主元或每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。
高斯列主消元数值分析实验报告
数值分析实验报告之高斯列主消元法一、实验目的：清楚高斯列主元消去法与高斯主元素消去法的区别，以及它提出的必要性；掌握高斯列主元消去法的原理及推导过程，会用其解简单的线性方程组。
二、实验内容：用高斯列主元消去法解方程组⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--643.5072.1000.2623.4712.3000.1000.3000.2001.0⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡321x x x =⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡000.3000.2000.1 三、实验原理：在采用高斯消去法解方程组时，小主元可能产生麻烦，即用其做除数，会导致其他元素数量级的严重增长和舍入误差的扩散，最终使得计算的解不可靠。
故应避免采用绝对值小的主元素。
在消元之前，选择一个绝对值最大的元素作为主元，用其做除数来进行消元，这样就具有较好的数值稳定性。
这就是选主元消去法。
下面详细说明列主元素消去法。
第一步：在Ax=b 即)1()1(b x A =的系数矩阵)1(A 的第一列元素中选择一个绝对值最大的元素，不妨设为)1(1l a 。
对调)1(1j a 和)1(lj a 及)1(1b 和)1(l b (j=1,2,……,n ，1≤l ≤n)。
以)1(1l a 作为新的)1(11a 进行消元（消去对调后的第2～n 个方程中的1x ）。
第k 步：（1≤k ≤n-1）设第k-1步消元过程完成，得到)()(k k b x A =，检查)(k A 中第k 列的后n-k+1个元素)(k kk a ，)(1k k k a +，…，)(k nk a ，从中选出绝对值最大者，不妨设是)(k pk a ，称它为第k 列主元素。
若p=k ，则取)(k kk a 做除数直接进行消元。
若p ≠k,则将第p 个方程与第k 个方程对调，使)(k pk a 成为新的)(k kk a ，然后以其作为除数进行消元，继续这一过程，直至得到等价的三角形方程组)()(n n b x A =，下一阶段的回代过程不变。
Gauss消去法实验报告
程序代码
程序代码说明
function x=gauss(A,b)
n=length(A);
a=[A,b];
for k=1:n-1
maxa=max(abs(a(k:n,k)));
if maxa==0
return;
end
for i=k:n
end
end
x=zeros(n,1);
x(n)=b(n)/a(n,n);
for k=n-1:-1:1
s=b(k);
for j=k+1:n
s=s-a(k,j)*x(j);
end
x(k)=s/a(k,k);
end
x
％调用M文件gauss. m
％形成增广矩阵
％计算乘子
％ 对k+1~n项进行消元
%增广矩阵第i行减去第k行的乘子倍目的是将该矩阵中的第k列中a(k,k)以下的元素全部消为零
二、实验原理：
由一般线性方程组在使用Gauss消去法求解时，从求解过程中可以清楚地看到，若 ，必须施以行交换的手续，才能使消去过程继续下去。有时既使 ，但其绝对值很小，由于舍入误差的影响，消去过程也会出现不稳定现象。因此，为使这种不稳定现象发生的可能性减至最小，在施行消去过程时每一步都要选主元素，即要寻找行 ，使
这样就完成了第1步消元。
回代过程：
在最后的一方程中解出 ，得：
再将 的值代入倒数第二个方程，解出 ，依次往上反推，即可求出方程组的解：
其通项为
三、实验内容与步骤：
１、实验内容：依照实验原理编写gauss消去法的程序。
２、实验步骤：首先，在电脑上安装matlab，然后，启动matlab，新建一个M文件。
2020年数值分析实验作业,gauss消去法的数值稳定性分析
作者：败转头作品编号44122544：GL568877444633106633215458 时间：2020.12.13实验3.1 Gauss 消去法的数值稳定性试验实验目的:观察和理解Gauss 消元过程中出现小主元(即)(k kka 很小)时引起的方程组解的数值不稳定性。
实验内容：求解方程组b Ax =，其中（1）⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⨯=11212592.1121-130.6-291.51314.59103.015-1A ，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=2178.4617.591b ； （2）⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----=2010151526990999999999.23107102A ,⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=15019000000000.582b . 实验要求：（1） 计算矩阵的条件数，判断系数矩阵是良态的还是病态的。
（2） 用Gauss 列主元消去法求得L 和U 及解向量421,R x x ∈.（3） 用不选主元的Gauss 消去法求得L ~和U ~及解向量421~,~R x x ∈.（4） 观察小主元并分析其对计算结果的影响.程序如下:计算矩阵条件数及Gauss列主元消去法: format longengA1=[0.3e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1];b1=[59.17;46.78;1;2];n=4;k2=cond(A1) %k2为矩阵的条件数；for k=1:n-1a=max(abs(A1(k:n,k)));[p,k]=find(A1==a);B=A1(k,:);c=b1(k);A1(k,:)=A1(p,:);b1(k)=b1(p);A1(p,:)=B;b1(p)=c;if A1(k,k)~=0A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k);A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n); elsebreakendendL1=tril(A1,0);for i=1:nL1(i,i)=1;endL=L1U=triu(A1,0)for j=1:n-1b1(j)=b1(j)/L(j,j);b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j);endb1(n)=b1(n)/L(n,n);for j=n:-1:2b1(j)=b1(j)/U(j,j);b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j);endb1(1)=b1(1)/U(1,1);x1=b1运行结果如下：K2=68.43；⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡---⨯=-14929.00202.00893.0011755.04724.00011079.26000118L ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=801.0000231.1835.2001314.5902592.11U 1x =[18.9882;3.3378;-34.747;-33.9865]不选主元的Gauss 消去法程序：clearformat longengA1=[0.3e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1]; b1=[59.17;46.78;1;2]; n=4;for k=1:n-1A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k);A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n); endL1=tril(A1,0); for i=1:nL1(i,i)=1; end L=L1U=triu(A1,0) for j=1:n-1b1(j)=b1(j)/L(j,j);b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j); endb1(n)=b1(n)/L(n,n); for j=n:-1:2b1(j)=b1(j)/U(j,j);b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j); endb1(1)=b1(1)/U(1,1); x1=b1程序运行结果如下：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⨯⨯=10189.010333.3011168.21033.3700110637.170001~151515L ⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡--⨯-⨯-⨯-⨯=-5.000816010637.171091.5210043.101314.59103.0~15151815U ]0;0;0005.1;6848.23[~1=x同理可得2A 对应的系数矩阵条件数及Gauss 列主元消去法求解结果： K2=8.994;⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-⨯=1333.04.0001104.0-3.0-0015.0000112-L ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=36667.300030.26005.155.2010710U ]0.1;;0.1;0.1;10444.0[152-⨯=-x不选主元的Gauss 消去法结果：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⨯=1400.0109999.0-001104998.2-5.00013.0-0001~1212L ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⨯⨯--=-3667.3000107495.5109987.14003.26100.1010710~121212U ]000145.1;99994.0;000.1;1045.1[~52-⨯-=-x实验4.5 三次样条插值函数的收敛性问题提出：作者：败转头作品编号44122544：GL568877444633106633215458时间：2020.12.13多项式插值不一定收敛的，即插值的节点多，效果不一定就好。
数值分析3 Gauss消去法
§１ Gauss 消去法它是以行初等变换为基础，整理而成的计算机算法。
由它改进和变形得到的高斯选主元消去法及三角分解法，仍然是目前计算机上常用的解低阶稠密矩阵的线性方程组的有效方法。
一、 引例：解方程组：⎪⎩⎪⎨⎧-=-+-=+-=+-4336425294321321321x x x x x x x x x用Matlab 先画出空间三个平面，并找出三平面的交点clear,clf,x=[-9 –9 9 9 -9];y=[-9 9 9 –9 -9];z=5-4*x+9*y;plot3(x,y,z),hold on,grid on z=3-2*x+4*y;plot3(x,y,z) z=-4+x-y;plot3(x,y,z)A=[4 –9 2;2 –4 6;-1 1 –3];b=[5 3 –4]';x=A\b; plot3(x(1),x(2),x(3),'r*')解：方程组的增广矩阵是[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-----=431136425294b A消元运算：[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--−−−→−⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----−−→−-----5.110005.055.00529475.25.225.105.055.0052942313125.025.14142r r r r r r b A回代求解：0.15,1.5/103-=-=x [],5.25.0/)15.0(55.0x 2=-⨯-=95.64/)15.025.295(x 1=⨯+⨯+=。
【注】Gauss 消元的基本思想：用初等行变换方法化一般的方程组为上三角方程组，然后求解；Gauss 消去法步骤：先消元计算，再回代求解二、Gauss 消去法设ｎ阶线性代数方程组为Ax=b (1)其中系数矩阵Ａ(n 阶非奇异矩阵)和右端列向量b 分别为A=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡nn n n n n a a a a a a a a a 212222111211 , b=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡n b b b 21 （一）上三角方程组的回代求解：⎪⎪⎩⎪⎪⎨⎧==++=+++n n nn n n n n b x u b x u x u b x u x u x u 2222211212111 ………………（2） 假设u ii ≠０（i=1,2,…，n ），则回代解得：⎪⎩⎪⎨⎧-==∑+=ii ni j j ij i innn n u x u b x u b x /)(/11,2,2,1 --=n n i ………………(3) （二）ＧＡＵＳＳ消元过程：(行变换)记方程组（1）为⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++)1()1(2)1(21)1(1)1(2)1(22)1(221)1(21)1(1)1(12)1(121)1(11nn nn n n n n n n b x a xa x ab x a x a x a b x a x a x a (4)对应的增广矩阵为[]⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=)1()1()1(2)1(1)1(2)1(2)1(22)1(21)1(1)1(1)1(12)1(11)1()1(,n nn n n nn b a a a b a a a b a a a bA…………………(5) 第一轮消元：若)1(11a ≠０，利用)1(11a 消下面的)1(1i a (i=2,3,…n), 消元因子：)1(11)1(11/a a l i i = (i=2,3,…n)[]⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=)2()2()2(2)2(2)2(2)2(22)1(1)1(1)1(12)1(11)2()2(00,n nn n nnb a a b a a b a a a bA………………(6) 其中)1(11)1()2(j i ij ij a l a a -= n j i ,,3,2, =)1(11)1()2(b l b b i i i -= n i ,,3,2 = (7)第二轮消元：若0)2(22≠a ，用第i 行减去第二行的)2(22)2(22/a a l i i =(i=3,4,…n)倍，得到[]⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=)3()3()3(3)3(3)3(3)3(33)2(2)2(2)2(23)2(22)1(1)1(1)1(13)1(12)1(11)3()3(00000,n nnn nn nb a a b a a b a a a b a a a a b A……………(8) 其中n j i b l b b a l a a i i ij i ij ij ,,4,3,,,)2(22)2()3()2(22)2()3( =-=-= …………(9) 已进行k-1轮消元后，得[]⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=)()()()()()()2(2)2(2)2(2)2(22)1(1)1(1)1(1)1(12)1(11)()(,k n k nnk nk k k k knk kk nkn k k K b a a b a a b a a a b a a a a b A……………(10) 第k 轮消元：若)(k kk a ≠0，用第i 行元素减去第k 行元素的)()(/k kk k ik ik a a l = (i=k+1,k+2,…,n)倍，得[]⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=++++++++++++++++)1()1()1(111)1(1)1(11)()()(1)()2(2)2(2)2(12)2(2)2(22)1(1)1(1)1(11)1(1)1(12)1(11)1()1(,k n k nnk nk k k k nk k k k k k kn knk kk k kknk kn k kk K b a a b a a b a a a b a a a a b a a a a a b A…(11) 其中⎪⎩⎪⎨⎧-=-=++,)()()1()()()1(k k ik k i k ik kjik k ij k ij b l b b a l a a ,,,1,n k j i += (12)最后，得[]⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=)()()2(2)2(2)2(22)1(1)1(1)1(12)1(11)()(,n n n nnnnn n b a b a a b a a a bA消元过程：[][][][])()()()()2()2()1()1(,,,,n n k k b A b A b A b A−→−−→−−→−−→−−→− （三）算法公式消元公式：对k=1,2,…,n-1若)(k kk a ≠0 )()(/k kk k ik ik a a l =⎪⎩⎪⎨⎧-=-=++,)()()1()()()1(k k ik k i k ik kj ik k ij k ij b l b b a l a a ,,,1,n k j i +=回代公式：)(n nn a ≠0⎪⎩⎪⎨⎧-==∑+=)(1)()()()(/)(/i ii n i j j i ij i i in nn n n n a x a b x a b x1,2,2,1 --=n n i三、GAUSS 消去法得以实现的前提(可行性分析)1、0a )()2(22(1)11≠n nn a a 或 2、,0,,0,,0,02122221112112221121111≠≠=≠≠A a a a aa a a a a D a a a a a kkk k k kk (13)这是因为：⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡→⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡)()2(2)2(22)1(1)1(12)1(11212222111211k kk n k kk k k k k a a a a a a a a a a a a a a a 行变换 k i a D n k i ii k ,,2,1,0,0,,,2,1)( =≠≠=则有若令【如】对称正定矩阵、严格对角占优矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡632351214 ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-103131021310四、Gauss 消去法的运算量消元过程乘法:[])1(31)()1(211-=-⨯+-∑-=n n k n k n n k 除法:)1(2)(11-=-∑-=n nk n n k 回代过程乘除法总运算量:),1(221+=+++n nnGauss 消去法的总运算量为:3)3(31)1(2)1(2)1(313232n n n n n n n n n n ≈-+=++-+- 对比克莱姆法则的总运算量为(n+1)!【注】对比：当n=10时，克莱姆法则(n+1)!=39916800，而Gauss 消去法，运算４３０次五、算法及例子步１，定义数组A,B,x ； 步２，输入A,B,Eps ；步３，对k=1到n-1进行消元：⎪⎩⎪⎨⎧++=-⇐>-⇐n k k j i a b a b b Eps a a a a a a kkk ik i i kk kk kj ik ij ij ,,2,1,/,,/ 步4，nn n n a b x /=；对i=n-1到1回代解ij ni j j iji i a x ab x /)(1∑+=-=步５，输出n i x i ,,2,1, =； 步６，结束。
数值分析实验(Gauss消元法)
end
%这里采用增广矩阵行变换的方式求解
c=n+1;
A(:,c)=b;
%%消元过程
for k=1:n-1
A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c);
end
%%回代结果
x=zeros(length(b),1);
x(n)=A(n,c)/A(n,n);
教师签名
年月日
重庆大学
学生实验报告
实验课程名称数值计算
开课实验室DS1421
学院2009年级土木专业班11班
学生姓名周波学号********
开课时间2010至2011学年第一学期
总成绩
教师签名
课程
名称
数值计算
实验项目
名称
Gauss消元法
实验项目类型
验证
演示
综合
设计
其他
指导
教师
何光辉
成绩
√
实验目的、实验原理
通过高斯消元法的实践环节，达到加深对高斯消元法理解和培养学生程序设计能力的目的。
在Matlab命令窗口中输入Gauss得出结果如下:
>> Gauss
x=
2.0000
1.0000
-1.0000
实验中遇到的问题及解决办法
消元过程要求a(ii)≠0,(i=1,2,…,n-1),回代过程则进一步要求a(nn)≠0。
Matlab求解线性方程的几种命令如下(方程组的一般形式可用矩阵和向量表示成: ,但运用下列方法的前提必须保证所求解的方程为恰定方程,即方程组存在唯一的一组解) :
运用求逆思想 : 或 ;
数值分析实验作业,gauss消去法的数值稳定性分析
实验3.1 Gauss 消去法的数值稳定性试验实验目的:观察和理解Gauss 消元过程中出现小主元(即)(k kka 很小)时引起的方程组解的数值不稳定性。
实验内容：求解方程组b Ax =，其中（1）⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⨯=11212592.1121-130.6-291.51314.59103.015-1A ，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=2178.4617.591b ； （2）⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----=2010151526990999999999.23107102A ,⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=15019000000000.582b .实验要求：（1） 计算矩阵的条件数，判断系数矩阵是良态的还是病态的。
（2） 用Gauss 列主元消去法求得L 和U 及解向量421,R x x ∈.（3） 用不选主元的Gauss 消去法求得L ~和U ~及解向量421~,~R x x ∈.（4） 观察小主元并分析其对计算结果的影响.程序如下:计算矩阵条件数及Gauss 列主元消去法:format longengA1=[0.3e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1]; b1=[59.17;46.78;1;2]; n=4;k2=cond(A1) %k2为矩阵的条件数；for k=1:n-1a=max(abs(A1(k:n,k))); [p,k]=find(A1==a); B=A1(k,:);c=b1(k);A1(k,:)=A1(p,:);b1(k)=b1(p); A1(p,:)=B;b1(p)=c; if A1(k,k)~=0A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k);A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n); else break end endL1=tril(A1,0); for i=1:n L1(i,i)=1; end L=L1U=triu(A1,0) for j=1:n-1b1(j)=b1(j)/L(j,j);b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j); endb1(n)=b1(n)/L(n,n); for j=n:-1:2b1(j)=b1(j)/U(j,j);b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j); endb1(1)=b1(1)/U(1,1); x1=b1运行结果如下： K2=68.43；⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡---⨯=-14929.00202.00893.0011755.04724.00011079.2600118L⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=801.0000231.1835.2001314.5902592.11U 1x =[18.9882;3.3378;-34.747;-33.9865] 不选主元的Gauss 消去法程序：clearformat longengA1=[0.3e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1]; b1=[59.17;46.78;1;2]; n=4;for k=1:n-1A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k);A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n); endL1=tril(A1,0); for i=1:n L1(i,i)=1; end L=L1U=triu(A1,0) for j=1:n-1b1(j)=b1(j)/L(j,j);b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j); endb1(n)=b1(n)/L(n,n); for j=n:-1:2b1(j)=b1(j)/U(j,j);b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j); endb1(1)=b1(1)/U(1,1); x1=b1程序运行结果如下：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⨯⨯=10189.010333.3011168.21033.3700110637.170001~151515L⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡--⨯-⨯-⨯-⨯=-5.000816010637.171091.5210043.101314.59103.0~15151815U ]0;0;0005.1;6848.23[~1=x同理可得2A 对应的系数矩阵条件数及Gauss 列主元消去法求解结果： K2=8.994;⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-⨯=1333.04.0001104.0-3.0-0015.0000112-L ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=36667.300030.26005.155.2010710U ]0.1;;0.1;0.1;10444.0[152-⨯=-x不选主元的Gauss 消去法结果：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⨯=1400.0109999.0-001104998.2-5.00013.0-0001~1212L ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⨯⨯--=-3667.3000107495.5109987.14003.26100.1010710~121212U ]000145.1;99994.0;000.1;1045.1[~52-⨯-=-x实验4.5 三次样条插值函数的收敛性问题提出：多项式插值不一定收敛的，即插值的节点多，效果不一定就好。
数值分析实验报告
实验五 解线性方程组的直接方法实验5.1 （主元的选取与算法的稳定性） 问题提出：Gauss 消去法是我们在线性代数中已经熟悉的。
但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss 消去法作为数值算法的稳定性呢？Gauss 消去法从理论算法到数值算法，其关键是主元的选择。
主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。
实验内容：考虑线性方程组n n n R b R A b Ax ∈∈=⨯,,编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss 消去过程。
实验要求：（1）取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816 b A ，则方程有解T x )1,,1,1(* =。
取n=10计算矩阵的条件数。
让程序自动选取主元，结果如何？（2）现选择程序中手动选取主元的功能。
每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。
若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。
（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。
（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。
重复上述实验，观察记录并分析实验结果。
思考题一：（Vadermonde 矩阵）设⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=∑∑∑∑====n i i n n i i ni i n i i n n n n n n nx x x x b x x x x x x x x x x x x A 002010022222121102001111 ，， 其中，n k k x k ,,1,0,1.01 =+=，（1）对n=2,5,8，计算A 的条件数；随n 增大，矩阵性态如何变化？（2）对n=5，解方程组Ax=b ；设A 的最后一个元素有扰动10-4，再求解Ax=b（3）计算（2）扰动相对误差与解的相对偏差，分析它们与条件数的关系。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析实验作业,g a u s s消去法的数值稳定性分析实验3.1 Gauss 消去法的数值稳定性试验实验目的:观察和理解Gauss 消元过程中出现小主元(即)(k kka 很小)时引起的方程组解的数值不稳定性。
实验内容：求解方程组b Ax =，其中（1）⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⨯=11212592.1121-130.6-291.51314.59103.015-1A ，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=2178.4617.591b ； （2）⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----=2010151526990999999999.23107102A ,⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=15019000000000.582b . 实验要求：（1） 计算矩阵的条件数，判断系数矩阵是良态的还是病态的。
（2） 用Gauss 列主元消去法求得L 和U 及解向量421,R x x ∈.（3） 用不选主元的Gauss 消去法求得L ~和U ~及解向量421~,~R x x ∈.（4） 观察小主元并分析其对计算结果的影响. 程序如下:计算矩阵条件数及Gauss 列主元消去法: format longengA1=[0.3e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1];b1=[59.17;46.78;1;2]; n=4;k2=cond(A1) %k2为矩阵的条件数；for k=1:n-1a=max(abs(A1(k:n,k))); [p,k]=find(A1==a); B=A1(k,:);c=b1(k);A1(k,:)=A1(p,:);b1(k)=b1(p); A1(p,:)=B;b1(p)=c; if A1(k,k)~=0A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k);A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n); else break end endL1=tril(A1,0); for i=1:n L1(i,i)=1; end L=L1U=triu(A1,0) for j=1:n-1b1(j)=b1(j)/L(j,j);b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j); endb1(n)=b1(n)/L(n,n); for j=n:-1:2b1(j)=b1(j)/U(j,j);b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j); endb1(1)=b1(1)/U(1,1); x1=b1运行结果如下： K2=68.43；⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡---⨯=-14929.00202.00893.0011755.04724.00011079.2600118L⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=801.0000231.1835.2001314.5902592.11U 1x =[18.9882;3.3378;-34.747;-33.9865] 不选主元的Gauss 消去法程序：clearformat longengA1=[0.3e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1]; b1=[59.17;46.78;1;2]; n=4;for k=1:n-1A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k);A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n); endL1=tril(A1,0); for i=1:n L1(i,i)=1; end L=L1U=triu(A1,0) for j=1:n-1b1(j)=b1(j)/L(j,j);b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j); endb1(n)=b1(n)/L(n,n); for j=n:-1:2b1(j)=b1(j)/U(j,j);b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j); endb1(1)=b1(1)/U(1,1); x1=b1程序运行结果如下：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⨯⨯=10189.010333.3011168.21033.3700110637.170001~151515L⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡--⨯-⨯-⨯-⨯=-5.000816010637.171091.5210043.101314.59103.0~15151815U ]0;0;0005.1;6848.23[~1=x同理可得2A 对应的系数矩阵条件数及Gauss 列主元消去法求解结果： K2=8.994;⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-⨯=1333.04.0001104.0-3.0-0015.0000112-L ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=36667.300030.26005.155.2010710U ]0.1;;0.1;0.1;10444.0[152-⨯=-x不选主元的Gauss 消去法结果：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⨯=1400.0109999.0-001104998.2-5.00013.0-0001~1212L ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⨯⨯--=-3667.3000107495.5109987.14003.26100.1010710~121212U ]000145.1;99994.0;000.1;1045.1[~52-⨯-=-x实验4.5 三次样条插值函数的收敛性问题提出：多项式插值不一定收敛的，即插值的节点多，效果不一定就好。
对三次样条插值函数又如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，也超过了本课程的内容。
通过本实验可以验证这一理论结果.实验内容：请按一定的规则分别选择等距或者非等距的插值节点，并不断增加插值节点的个数。
考虑实验4.4中的函数或者选择其他感兴趣的函数，可用Matlab的函数“spline”作此函数的三次样条插值函数。
实验要求：（1）随着节点个数的增加，比较被逼近函数和三次样条差值函数的误差变化情况。
分析所得结果并与拉格朗日插值多项式比较。
（2）三次样条插值函数的思想最早产生于工业部门。
作为工业迎合用的例子，考虑如下例子：某汽车制造商根据三次样条差值函数设计车门曲线，其中一段的数据如下：(3)计算实验4.4的样条插值.程序如下：format shortx1=-1:0.5:1;y1=1./(1+25.*x1.*x1);x2=-1:0.25:1;y2=1./(1+25.*x2.*x2);x3=-1:0.1:1;y3=1./(1+25.*x3.*x3);x4=[-1,-0.82,-0.6,-0.53,-0.34,-0.2,0,0.04,0.2,0.25,0.5,0.8,1];y4=1./(1+25.*x4.*x4);xx=-1:0.01:1;yy1=spline(x1,y1,xx);yy2=spline(x2,y2,xx);yy3=spline(x3,y3,xx);yy4=spline(x4,y4,xx);hold onfplot('1./(1+25.*x.*x)',[-1,1],'m')plot(xx,yy1,'g')plot(xx,yy2,'b')plot(xx,yy3,'k')plot(xx,yy4,'r')hold off %比较被逼近函数与三次样条插值函数图像，直观表现不同插值节点处误差的变化xx=-1:0.2:1;y=1./(1+25.*xx.*xx)%函数在相应节点处的真实值；yy1=spline(x1,y1,xx)y1la=lagrange(x1,y1,xx)yy2=spline(x2,y2,xx)y2la=lagrange(x2,y2,xx)yy3=spline(x3,y3,xx)y3la=lagrange(x3,y3,xx)yy4=spline(x4,y4,xx)y4la=lagrange(x4,y4,xx)其中lagrange函数对应的m文件为：function y=lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;p=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;end程序运行结果如下：插值结果比较：y 0.0385 0.0588 0.1 0.2 0.5 1 0.5 0.2 0.1 0.0588 0.0385 yy1 0.0385 -0.252 -0.0623 0.3687 0.8024 1 0.8024 0.3687 -0.0623 -0.252 0.0385 y1la 0.0385 -0.3793 -0.1101 0.4005 0.8342 1 0.8342 0.4005 -0.1101 -0.3793 0.0385 yy2 0.0385 0.0551 0.1085 0.1765 0.5342 1 0.5342 0.1765 0.1085 0.0551 0.0385 y2la 0.0385 -0.2587 0.3049 0.0726 0.5636 1 0.5636 0.0726 0.3049 -0.2587 0.0385 yy3 0.0385 0.0588 0.1 0.2 0.5 1 0.5 0.2 0.1 0.0588 0.0385 y3la 0.0385 0.0588 0.1 0.2 0.5 1 0.5 0.2 0.1 0.0588 0.0385 yy4 0.0385 0.0587 0.1 0.2016 0.5 1 0.5 0.1989 0.0993 0.0588 0.0385 y4la 0.0385 0.4312 0.1 0.2461 0.5 1 0.5 0.264 0.2387 0.0588 0.0385 车门曲线求解程序:x=0:10;y=[0,0.79,1.53,2.19,2.71,3.03,3.27,2.89,3.06,3.19,3.29];dx10=0.2;S=csfit(x,y,dx0,dx10)其中csfit函数的m文件为:function S=csfit(X,Y,dx0,dxn)%Clamped Cubic Spline%Input -X is the 1xn abscissa vector% -Y is the 1xn ordinate vector% -dx0=S'(x0) first derivative boundary condition% -dxn=S'(xn) first derivative boundary condition%Output-S: rows of S are the coefficients, in descending order, for the % cubic interpolantsN=length(X)-1;H=diff(X);D=diff(Y)./H;A=H(2:N-1);B=2*(H(1:N-1)+H(2:N));C=H(2:N);U=6*diff(D);%Clamped spline endpoint constraintsB(1)=B(1)-H(1)/2;U(1)=U(1)-3*(D(1)-dx0);B(N-1)=B(N-1)-H(N)/2;U(N-1)=U(N-1)-3*(dxn-D(N));for k=2:N-1temp=A(k-1)/B(k-1);B(k)=B(k)-temp*C(k-1);U(k)=U(k)-temp*U(k-1);endM(N)=U(N-1)/B(N-1);for k=N-2:-1:1M(k+1)=(U(k)-C(k)*M(k+2))/B(k);endM(1)=3*(D(1)-dx0)/H(1)-M(2)/2;M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;for k=0:N-1S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));S(k+1,2)=M(k+1)/2;S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6;S(k+1,4)=Y(k+1);endend程序运行结果为:S=(3)实验4.4中的样条插值见上表中的y值.。