矩阵与数值分析(完整版)实用资料

矩阵与数值分析（完整版）实用资料(可以直接使用，可编辑完整版实用资料，欢迎下载）

矩阵与数值分析矩阵与数值分析（完整版）实用资料(可以直接使用，可编辑完整版实用资料，欢迎下载）学院专业班级学号姓名电子信息与电气工程学部生物医学工程刘江涛1：考虑计算给定向量的范数；输入向量x=(x1,x2,,xn)T，输出x,x2,x∞，请编制一个通用程序，并用你编制的程序计算如下向量的范数：1⎫⎛11Tx=1,,,,⎪,y=(1,2,,n)n⎭⎝23对n=10,100,1000甚至更大的n计算其范数，你会发现什么结果？你能否修改你的程序使得计算结果相对精确呢？通用求范数程序：functionNORM(x)y1=sum(abs(x));y2=(sum(x.^2))^(1/2);y3=max(abs(x));fprintf('1-范数=%g；2-范数=%g；inf-范数=%g\n',y1,y2,y3);例题的运行程序：functionxianglaing(n)x=[];y=[];fori=1:nx(i)=1/i;y(i)=i;enddisp('x的范数：');NORM(x');disp('')disp('y的范数：');NORM(y');运行结果如下表：T根据上述的两个表的运行结果，我们可以得知无论n的值如何变化，对于x∞=1恒成立；y∞=n恒成立，其1-范数与2-范数随着n的增大而增大，但是其变化越来越小，这是因为计算在进行数值计算时有误差存在，对于表达式(1)当n很大时1却很n小，会出现“大数吃小数的现象”；修改方案：当n很大时我们避免用n做除数，因为当n非常大时1→0成立；所以在求解其范数时我们从小数开始相加，无穷个非常n小的数值相加也可能是个很大的数，从而可以避免两个数相加时出现“大数吃小数”的现象；2：考虑y=f(x)=ln(1+x)，其中定义f(0)=1，此时f(x)是连续函数，用此公x式计算当x∈[-10-15,10-15]时的函数值，画出图像。另一方面，考虑下面算法：d=1+x;ifd=1theny=1elsey=lnd/(d-1)endif用此算法计算x∈[-10-15,10-15]时的函数值，画出图像，比较一下发生了什么？程序：x=-10^(-15):10^(-20):10^(-15);if(x==0)f=1;elsef=log(1+x)/x;endfigure(1)plot(x,f);d=1+x;ifd==1y=1;elsey=log(d)/(d-1);endfigure(2)Plot(x,y);有图可知，直接用公式f(x)=ln(1+x)计算x∈[-10-15,10-15]的函数值时，除了在xx=0出的值为1，其他的值都是无限趋近于1；而利用算法二算出的结果全为1；出现这这情况的原因是x的取值非常接近于0，在用公式d=1+x求d得过程中出现了大数吃小数的情况，所以在用计算机计算时d=1恒成立，从而使y=1恒成立；3：首先编写一个利用秦九韶算法计算一个多项式在定点的函数值的通用程序，你的程序包括输入多项式的系数以及定点，输出函数值，利用你编写的程序计算f(x)=(x-2)9=x9-18x8+144x7-672x6+2021x5-4032x4+5376x3-4608x2+2304x-512在x=2邻域附近的值，画出p(x)在x∈[1.95,20.5]上地图像。秦九韶算法的通用程序：%A为多项式的以升幂排列的系数，x为初始值functionp=qinjiushao(A,x)a=A;[~,n]=size(a);n=n-1;S=[];S(n+1)=a(n+1);fork=n:-1:1S(k)=x.*S(k+1)+a(k);endp=S(1);利用上述程序计算p(x)在x=2邻域附近的值具体见下表：当x∈[1.95,20.5]时，p(x)的图像如下：画图程序如下：functionhuatu(A,x)[~,n]=size(x);fori=1:ny(i)=qinjiushao(A,x(i));endplot(x,y);程序运行如下：>>x=1.95:0.01:20.5;>>A=[-5122304-46085376-40322021-672144-181];>>huatu(A,x)11p(x)x4：编制计算机给定矩阵A的LU分解和PLU分解的通用程序，然后利用你编写的程序完成下面两个计算任务：考虑⎡1⎢-1⎢A=⎢⎢⎢-1⎢⎣-101⎤⎥⎥01⎥∈Rn⨯n⎥11⎥-11⎥⎦0-1-1自己取定x∈Rn，并计算b=Ax。然后用你编制的不选主元的Gauss消去法求解ˆ。对n从5到30估计计算解的精度。该方程组，记你计算出的解为x（2）对n从5到30计算出其逆矩阵。LU分解的通用程序：functionLU(A)[m,n]=size(A);L=zeros(m,n);U=zeros(m,n);forj=1:nU(1,j)=A(1,j);L(j,j)=1;endforj=2:mL(j,1)=A(j,1)/U(1,1);endfori=2:nforj=i:nsum=0;fork=1:i-1sum=sum+L(i,k)*U(k,j);endU(i,j)=A(i,j)-sum;endforj=i+1:nsum=0;fork=1:i-1sum=sum+L(j,k)*U(k,i);endL(j,i)=(A(j,i)-sum)/U(i,i);endenddisp('L=');disp(L);disp('U=');disp(U);PLU分解的通用程序：functionPLU(A)[~,n]=size(A);Ip=1:n;fork=1:n-1[~,r]=max(abs(A(k:n,k)));r=r+(k-1);ifr>kA([k,r],:)=A([r,k],:);Ip([k,r])=Ip([r,k]);endforp=k+1:nmu=A(p,k)/A(k,k);A(p,k)=mu;A(p,k+1:n)=A(p,k+1:n)-mu*A(k,k+1:n);endendp=eye(n,n);P=zeros(n,n);fori=1:nP(i,1:n)=p(Ip(i),1:n);endL=tril(A,-1)+eye(n);U=triu(A);disp('P=')disp(P);disp('L=')disp(L);disp('U=');disp(U);我选取b=[123n]Tn∈[5,30],n∈N*，则我们可以计算出其精确解2n-1-12n-2-1-1x=[-n-1-n-22222n-1T]2n-1n∈[5,30]n∈N*实现求解的程序如下：functionINV(n)A=ones(n);fori=2:nforj=1:i-1A(i,j)=-1*A(i,j);endforj=i:n-1A(i-1,j)=0;endendfori=1:nb(i)=i;X(i)=-(2^(n-i)-1)/(2^(n-i));endX(n)=(2^n-1)/(2^(n-1));[L,U]=LU(A);x1=inv(U)*inv(L)*b'disp(x1);我们利用上述的程序计算出的结果如下表：利用我们求出的精确解与用程序求出的近似解求其误差，再利用matlab编程实现时其误差为零。（2）对于题目中的A的求逆的通用程序：functionINV(n)A=ones(n);fori=2:nforj=1:i-1A(i,j)=-1*A(i,j);endforj=i:n-1A(i-1,j)=0;endendB=inv(A);disp('A的逆为：');disp(B);n=5时：-0.2500-0.1250-0.0625-0.0625⎤⎡0.5000⎢⎥00.5000-0.2500-0.1250-0.1250⎢⎥-1⎥A=⎢000.5000-0.2500-0.2500⎢⎥0000.5000-0.5000⎢⎥⎢⎥0.25000.12500.06250.0625⎣0.5000⎦n=6时：-0.2500-0.1250-0.0625-0.0313-0.0313⎤⎡0.5000⎢0⎥0.5000-0.2500-0.1250-0.0625-0.0625⎢⎥⎢0⎥00.5000-0.2500-0.1250-0.1250A-1=⎢⎥000.5000-0.2500-0.2500⎢0⎥⎢0⎥0000.5000-0.5000⎢⎥0.50000.25000.12500.06250.03130.0313⎢⎥⎣⎦n=7时：-0.2500-0.1250-0.0625-0.0313-0.0156-0.0156⎤⎡0.5000⎢0⎥0.5000-0.2500-0.1250-0.0625-0.0313-0.0313⎢⎥⎢0⎥00.5000-0.2500-0.1250-0.0625-0.0625⎢⎥A-1=⎢0000.5000-0.2500-0.1250-0.1250⎥⎢00000.5000-0.2500-0.2500⎥⎢⎥000000.5000-0.5000⎢⎥⎢0.50000.25000.12500.06250.03130.01560.0156⎥⎣⎦n=8时：-0.2500-0.1250-0.0625-0.0313-0.0156-0.0078-0.0078⎤⎡0.5000⎢0⎥0.5000-0.2500-0.1250-0.0625-0.0313-0.0156-0.0156⎢⎥⎢0⎥00.5000-0.2500-0.1250-0.0625-0.0313-0.0313⎢⎥0000.5000-0.2500-0.1250-0.0625-0.0625⎥A-1=⎢⎢00000.5000-0.2500-0.1250-0.1250⎥⎢⎥000000.5000-0.2500-0.2500⎢⎥⎢0000000.5000-0.5000⎥⎢⎥0.25000.12500.06250.03130.01560.00780.0078⎥⎢⎣0.5000⎦.......在此不再一一列举，都可以用上述程序算出；5：编制计算对称正定阵的Cholesky分解的通用程序，并利用你编制的程序计算Ax=b，其中A=(aij)∈Rn⨯n,aij=1，b可以有你自己取定，对n从10到20验i+j-1证程序的可靠性。Cholesky分解求L的通用程序：%LT代表L的转置function[L,LT]=Cholesky(A)[n,m]=size(A);L=zeros(n,n);forj=1:msum=0;fork=1:j-1sum=sum+(L(j,k))^2;endL(j,j)=(A(j,j)-sum)^(1/2);fori=j+1:nsum=0;fork=1:j-1sum=sum+L(i,k)*L(j,k);endL(i,j)=(A(i,j)-sum)/L(j,j);endendL=L;LT=L';求解Axb的方程解的程序如下：functioncholesky_qiu_jie(n,b)A=zeros(n);fori=1:nforj=1:nA(i,j)=1/(i+j-1);endend[L,LT]=Cholesky(A);Y=inv(LT)*inv(L)*b;disp(Y');>>n=10;>>b=[1000000000]';>>cholesky_qiu_jie(n,b)1.0e+06*Columns1through5Columns6through10-6.3052258749891819.607833253959356-8.7498890229612894.374900125373453-0.923581794866949>>n=11;>>b=[10000000000]';>>cholesky_qiu_jie(n,b)1.0e+07*Columns1through5Columns6through10Column110.385801129112326n=12;>>b=[111111111111]';>>cholesky_qiu_jie(n,b)1.0e+08*Columns1through5Columns6through10Columns11through12n=13;>>b=[1111111111111]';>>cholesky_qiu_jie(n,b)1.0e+10*Columns1through5Columns6through10-0.0955644689000000.338540141400000-0.795937993200000Columns11through13>>n=14;>>b=[11111111111110]';>>cholesky_qiu_jie(n,b)1.0e+15*Columns1through5Columns6through10Columns11through14-0.862636505012200>>n=15;>>b=[111111111111101]';>>cholesky_qiu_jie(n,b)1.0e+16*Columns1through5Columns6through10Columns11through15>>n=16;>>b=[1111111111111010]';>>cholesky_qiu_jie(n,b)1.0e+16*Columns1through5Columns6through10Columns11through15Column16>>n=17;>>b=[11111111111110101]';>>cholesky_qiu_jie(n,b)1.0e+17*Columns1through5Columns6through10-0.034029369522343Columns11through15Columns16through17>>n=18;>>b=[111111111111101012]';>>cholesky_qiu_jie(n,b)1.0e+16*Columns1through5Columns6through10Columns11through15Columns16through18-3.4602807804032691.887745959501761-0.436632060441746>>n=19;>>b=[1111111111111010120]';>>cholesky_qiu_jie(n,b)1.0e+17*Columns1through5Columns6through10Columns11through15Columns16through191.944997967101147-0.384572610629205n=20;b=[11111111111110101200]';>>cholesky_qiu_jie(n,b)1.0e+17*Columns1through5Columns6through10Columns11through15Columns16through203.295476378249488A-1利用cholesky分解求出的解的精确度高于直接，因为当n逐渐增大的过程中，越来越接近奇异矩阵，使得计算结果的误差增大，而使用cholesky分解可以避免这种现象的产生，是计算结果更加精确。6：（1）编制程序House(x)，其作用是对输入的向量x，输出单位向量u使得(I-2uuT)x=x2e1。编制Householder变换阵H=I-2uuT∈Rn⨯n乘以A∈Rn⨯m的程序HA，注意，你的程序并不显式的计算出H。考虑矩阵⎛1-1A=-2-0⎝232224⎫⎪2⎪eπ⎪⎪-37⎪75/2⎪⎭3用你编制的程序计算H使得HA的第一列为αe1，并将HA的结果显示出来。编制House(x)，其作用是对输入的向量x,输出单位向量U使得(I-2uuT)x=x2e1：House(x)通用程序：functionHouse(x)[n,~]=size(x);e1=zeros(n,1);e1(1)=1;w=x-norm(x)*e1;U=w/sqrt((w'*w));disp('U=');disp(U);编制Householder变换阵H=I-2uuT∈Rn⨯n乘以A∈Rn⨯n的程序HA，注意，你的程序并不显示的计算出H。Householder变换阵H通用程序：functionHouseholder(A)[n,~]=size(A);e1=zeros(n,1);e1(1)=1;x=A(:,1);w=x-norm(x)*e1;U=w/sqrt((w'*w));H=eye(n)-2*(U*U');HA=H*A;disp('HA=');disp(HA);考虑矩阵⎛1-1A=-2-0⎝24⎫⎪32⎪2eπ⎪⎪2-37⎪272⎪⎭3用你编制的程序计算H使得HA的第一列为αe1的形式，并将HA的结果显示：HA的结果显示的通用程序：functionHouseholder(A)[n,~]=size(A);e1=zeros(n,1);e1(1)=1;x=A(:,1);w=x-norm(x)*e1;U=w/sqrt((w'*w));H=eye(n)-2*(U*U');HA=H*A;disp('HA=');disp(HA);程序运行结果：>>A=[1234;-13sqrt(2)sqrt(3);-22exp(1)pi;-sqrt(10)2-37;0275/2];>>Householder(A)HA=4.0000-2.83111.4090-6.53780.00001.38960.8839-1.78050.0000-1.22081.6576-3.88360.0000-3.0925-4.6770-4.107802.00007.00002.50007：用jacobi和Gauss-Scidel迭代求解下面的方程组，输出每一步的误差xk-x*；⎧5x1-x2-3x3=-2⎪⎨-x1+2x2+4x3=1⎪-3x+4x+15x=10123⎩Jacobi迭代通用程序：functionjacobi(A,b,x0,ep)ifnargin==3ep=1.0e-6;elseifnargin<3errorreturnendendD=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=D\(L+U);f=D\b;y=B*x0+f;n=1;whilenorm(y-x0)>=epfprintf('%d步误差\n',n)Error=y-x0;disp(Error)x0=y;y=B*x0+f;n=n+1;enddisp('近似解')disp(y)在命令窗口输入以下命令就可以输出每一步的误差及最优近似解：>>A=[5-1-3;-124;-3415];>>b=[-2110]';>>x0=[000]';>>ep=10e-8;>>jacobi(A,b,x0,ep)近似解Gauss-Scidel迭代法通用程序：functionGaussseidel(A,b,x0,ep)ifnargin==3ep=1.0e-6;elseifnargin<3errorreturnendendD=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=(D-L)\U;f=(D-L)\b;y=B*x0+f;n=1;whilenorm(y-x0)>=epfprintf('%d步误差\n',n)Error=y-x0;disp(Error)x0=y;y=B*x0+f;n=n+1;enddisp('x的近似最优解：')disp(y);在命令窗口输入以下命令就可以输出每一步的误差及最优近似解：>>A=[5-1-3;-124;-3415];>>b=[-2110]';>>x0=[000]';>>ep=10e-8;>>Gaussseidel(A,b,x0,ep)x的近似最优解：-1.8032785851349118：取不同的初值用Newton迭代法以及弦截法求解x3+2x2+10x-100=0的实根，列表或者画图说明收敛性；Newton迭代法通用程序：%使用说明f为符号表达式，x0为初始值functionNewton(f,x0)symsx;x=x0;x1=x0-eval(f/diff(f));while(norm(x1-x0)>0.000001)x=x1;xk=x1-eval(f/diff(f));x0=x1;x1=xk;disp(xk);enddisp('xk=')disp(xk)弦截法通用程序：%使用说明f为符号表达式，x0,x1为初始值functionXuanjiefa(f,x0,x1)symsx;x=x0;f1=eval(f);x=x1;f2=eval(f);x2=x1-(f2/(f2-f1))*(x1-x0);while(norm(x2-x1)>0.000001)x=x1;f1=eval(f);x=x2;f2=eval(f);xk=x2-(f2/(f2-f1))*(x2-x1);x1=x2;x2=xk;enddisp('xk=')disp(xk)Newton法：弦截法：从上述两个表的结果我们可以对照得出，在不同的初始值的条件下运行程序都具有收敛性，但是不同的初始值有不同的收敛速度；9：用二分法求解方程excosx+2=0在区间[0,4π]上所有根。functioner_fen_fa(f,a,b)symsx;e=0.0000001;while((b-a)>e)m=(b+a)/2;x=a;f1=eval(f);x=m;f2=eval(f);iff1*f2>=0a=m;elseb=m;endendc=(a+m)/2;disp(c);根据函数的特点将其定义域分为4等分分别为[0,π]，[π,2π]，[2π,3π]，[3π,4π]分别运行程序：>>symsx>>f=exp(x)*cos(x)+2;>>a=0;>>b=pi;>>er_fen_fa(f,a,b近似零点：1.8807>>symsx>>f=exp(x)*cos(x)+2;a=pi;>>b=2*pi;>>er_fen_fa(f,a,b)近似零点：4.6941>>symsx>>f=exp(x)*cos(x)+2;>>a=2*pi;>>b=3*pi;>>er_fen_fa(f,a,b)近似零点：7.8548>>symsx>>f=exp(x)*cos(x)+2;>>a=3*pi;>>b=4*pi;>>er_fen_fa(f,a,b)近似零点：10.9955;x2=4.6941;x2=7.8548;x4=10.9955由上述程序可知其全部解为x1=1.880710：考虑函数f(x)=sin(πx),x∈[0,1]。用等距节点作f(x)的Newton插值，求出插值多项式以及f(x)的图像，观察收敛性。Newton插值求插值多项式的通用程序：functionP=Nweton_xhazhi(X,f)symsx;p=0;[~,n]=size(X);fori=1:n-1F=0;w=1;sum=1;forl=1:i+1;w=w*(x-X(l));endforj=1:i+1x=X(j);F=F+f(j)/(eval(diff(w)));endfort=1:isymsx;sum=sum*(x-X(t));endp=p+F*sum;endP=p+f(1);Disp(‘P=’);用等距节点作f(x)的Newton插值多项式的通用程序：X=0:0.1:1;f=sin(pi*X);O=Nweton_xhazhi(X,f);Y=simple(O);disp(Y);plot(X,f)P=画图程序：X=0:0.1:1;x1=0:0.0001:1;f1=sin(x1*pi);symsxx0=0:0.001:1;f=sin(pi*X);O=Nweton_xhazhi(X,f);x=x0;Y=eval(O);plot(x0,Y,'r');holdonplot(x1,f1);由上图可知：红色图形是利用插值多项式画出的，蓝色图像是利用原函数画出的，将两者对比我们可以得出结论：该插值多项式具有收敛性11：对函数f(x)=1,x∈[-5,5]，取不同的节点数n,用等距节点作Lagrange1+x2插值，观察Runge现象。Lagrange插值通用程序：%X为节点，y为节点对应的函数值，x0为插值节点functionLagrange(X,y,x0)symsx;[~,n]=size(X);sum=0;fori=1:nsum=sum+y(i)*L(X,i);endx=x0;p=eval(sum);plot(x0,p,'r-')holdonplot(X,y,'b-');functionl=L(X,k)symsx;sum1=1;sum2=1;[~,n]=size(X);if(k==1)fori=k+1:nsum1=sum1*(x-X(i));sum2=sum2*(X(k)-X(i));l=sum1/sum2;endelsefori=1:k-1sum1=sum1*(x-X(i));sum2=sum2*(X(k)-X(i));endfori=k+1:nsum1=sum1*(x-X(i));sum2=sum2*(X(k)-X(i));endl=sum1/sum2;end取不同的节点数n,用等距节点作Lagrange插值，观察Runge现象:当n=10;X=-5:0.1:5;f=1./(1+X.^2);x0=-5:5:5;Lagrange(X,f,x0)x1=-5:0.0001:5;y1=1./(1+x1.^2);plot(x1,y1)1.41.20.80.60.40.2当n=20;X=-5:0.1:5;f=1./(1+X.^2);x0=-5:0.5:5;Lagrange(X,f,x0)x1=-5:0.0001:5;y1=1./(1+x1.^2);plot(x1,y1)当n=50;X=-5:0.1:5;f=1./(1+X.^2);x0=-5:0.2:5;Lagrange(X,f,x0)x1=-5:0.0001:5;y1=1./(1+x1.^2);plot(x1,y1)1.41.210.80.60.40.2n=100;X=-5:0.1:5;f=1./(1+X.^2);x0=-5:0.1:5;Lagrange(X,f,x0)x1=-5:0.0001:5;y1=1./(1+x1.^2);plot(x1,y1)-5-4-3-2-1012345由以上几个图形我们可以得知，插值节点越多时，插值多项式的收敛性越好；12：令f(x)=e3xcos(πx)，考虑积分⎰2π0f(x)dx。区间分为50，100，200，500，1000等，分别用复合梯形以及复合Simpson积分公式计算积分值，将数值积分的结果与精确值比较，列表说明误差的收敛性。复合梯形公式通用程序：%f为被积函数，n为积分区间等分的数目，a,b分别为积分下上限functiontixing(f,n,a,b)symsx;sum=0;xk=a:(b-a)/n:b;fori=2:n-1x=xk(i);sum=sum+eval(f);symsx;endx=a;f1=eval(f);x=b;f2=eval(f);T=(b-a)/(2*n)*(f1+f2+2*sum);disp(T)复合sipson公式：functionSimpson(f,n,a,b)symsx;sum=0;sum1=0;xk=a:(b-a)/n:b;fori=1:n-1x=(xk(i)+xk(i+1))/2;sum=sum+eval(f);symsx;endfori=2:n-1x=xk(i);sum1=sum1+eval(f);symsx;endx=a;f1=eval(f);x=b;f2=eval(f);T=(b-a)/(6*n)*(f1+f2+2*sum1+4*sum);disp(T)利用matlab内置函数int()求出该积分式的真是值为：>>symsx>>int(exp(3*x)*cos(pi*x),0,2*pi)ans=3.5232e+07利用符合simpson公式计算如下表：利用复合梯形公式计算如下表：由上述两个表格的值我们可以得出，其误差随着n的增大近似值与真实值之间的误差逐渐减少，说明了误差具有收敛性。13：分别用2点,3点以及5点的Gauss型积分公式计算如下定积分：(1)⎰1x2-x2-1dx(2)⎰2πsinxdxx解：利用Gauss型求积公式得出：2点Gauss积分公式求解时，n=1;查找上表得x0,x1,A0,A1，把他们都代入，可得⎰1x2-x2-1dx≈(0.5773502692)2-(0.5773502692)2+(-0.5773502692)2-(-0.5773502692)2=0.81653点Gauss积分公式求解时，n=2；查找上表可得x0,x1,x2,A0,A1,A2，把他们都代入，可得：⎰1x2-x2-1dx≈0.5555555556⨯(0.7745966692)2-(0.7745966692)222+0.5555555556⨯(-0.7745966692)-(-0.7745966692)=1.05415点Gauss积分公式求解时，n=4;查找上表可得x0,x1,x2,x3,x4,A0,A1,A2,A3,A4，把他们都代入，可得：⎰1x22-1-x0.53846931012dx≈2⨯0.2369268851⨯=1.24952+2⨯0.4786286705⨯-0.53846931012⎰2π02πsinx1sin(π(t+1))sinx2π-02π+0dx，令x=t+=πt+π则⎰dx=⎰，于是利用0-1x22xt+1上表计算Gauss积分2点Gauss积分公式求解时，n=1;查找上表得x0,x1,A0,A1，把他们都代入，可得1sin(π(t+1))sinxsin(π(1+0.5773502692))dx==+⎰0x⎰-1t+11+0.5773502692sin(π(1-0.5773502692))=1.68121-0.57735026922π3点Gauss积分公式求解时，n=2；查找上表可得x0,x1,x2,A0,A1,A2，把他们都代入，可得：1sin(π(t+1))sinxsin(π(1+0.7745966692))dx==0.5555555556⨯⎰0x⎰-1t+11+0.7745966692sin(π(1-0.7745966692))+0.5555555556⨯=1.39951-0..77459666922π5点Gauss积分公式求解时，n=4;查找上表可得x0,x1,x2,x3,x4,A0,A1,A2,A3,A4，把他们都代入，可得：1sin(π(t+1))sinx51⨯⎰0xdx=⎰-1t+1dt=0.2369268859))sin(π(1-0.9061798459))⎫⎛sin(π(1+0.90617984+⎪+1+0.90617984591-0.9061798459⎝⎭2π01))sin(π(1-0.5384693101))⎫⎛sin(π(1+0.538469310.4786286705⨯+⎪1+0.53846931011-0.5384693101⎝⎭=1.418114：考虑微分方程初值问题：1⎧dx=(tx-x2)⎪2⎨dt(t+1)⎪x(0)=2⎩分别用Euler法，改进的Euler法，Runge-Kutta法求解该方程。分别去步长为0.1，0.01，0.001，计算到x(1),画图说明结果。Euler法：%f为被积表达式，h为步长，ab为t的取值范围，x0为初始值functionEuler(f,h,a,b,x0)symstxN=abs(b-a)/h;x=x0;Y=zeros(1,N+1);T=zeros(1,N+1);Y(1)=x0;forn=0:Nt=a+n*h;T(n+1)=t;if(n<N)un=x+h*eval(f);x=un;Y(n+2)=un;endenddisp(un);plot(T,Y,'*-')N=0.1;x(1)=1.072393608200700改进的Euler法：%f为被积表达式，h为步长，ab为t的取值范围，x0为初始值functionGai_Jin_Euler(f,h,a,b,x0)symstxN=abs(b-a)/h;x=x0;Y=zeros(1,N+1);T=zeros(1,N+1);Y(1)=x0;forn=0:Nt=a+n*h;T(n+1)=t;f1=eval(f);t=a+(n+1)*h;un1=x+h*f1;x=un1;f2=eval(f);ifn<Nx=x0;un=x+h/2*(f1+f2);Y(n+2)=un;x0=un;endx=un;enddisp(un);plot(T,Y,'*-');Runge-Kutta法：第33卷第3期岩土力学Vol.33No.32021年3月RockandSoilMechanicsMar.2021收稿日期:2021-09-20文章编号:1000-7598(202103-0906-07旋喷群桩复合地基承载特性的数值分析安关峰,张洪彬,刘添俊(广州市市政集团,广州510060摘要:旋喷桩加固软土地基在各种地基处理工程中得到了广泛应用。对旋喷桩的研究多数集中在其施工工艺的改进上,或者针对单桩的承载特性进行研究,而对旋喷群桩的承载特性则研究不多。根据工程实际情况,采用基于MIDAS-GTS的三维有限元分析技术,通过改变旋喷群桩的布置方式、桩弹性模量、桩长、桩径、桩距等设计参数及桩-土接触面等参数对旋喷群桩复合地基承载特性的影响进行了研究。研究表明:旋喷桩加固软土地基主要减小了地表至桩底深度范围内土体的竖向沉降,对桩底下方的土体沉降基本无影响;提高旋喷桩桩径及材料强度会提高复合地基承载能力;不同旋喷桩布置方式、桩-土之间是否设置Goodman接触面单元对地基承载能力基本无影响。关键词:旋喷桩;复合地基;承载特性;三维数值分析中图分类号:TU473.1文献标识码:ANumericalanalysisofbearingcharacteristicsofcompositesubgradereinforcedbychemicalchurningpilegroupsANGuan-feng,ZHANGHong-bin,LIUTian-jun(GuangzhouMunicipalEngineeringGroup,Guangzhou510060,ChinaAbstract:Themethodofsoftsoilsubgradereinforcementwithchemicalchurningpileisusedmoreandmoreinprojectsofsubgradetreatment.Atpresent,theresearchesaboutchemicalchurningpileusuallyfocusonimprovementinconstructiontechniquesorbearingcharacteristicsofsinglepile.Butresearchaboutbearingcharacteristicsofchemicalchurningpilegroupsisnotmuch.BasedonMIDAS-GTSthree-dimensionalfiniteelementanalysis,theinfluenceofdesignparametersonbearingcharacteristicsofchemicalchurningpilegroupsisstudied.Thesedesignparametersofchemicalchurningpileincludelayout,elasticmodulus,length,diameterofpiles,anddistancebetweenpiles,theparametersofinterfacebetweenpilesandsoil.Theresultsshowthattheverticalsettlementofsoilwithintherangebetweensurfaceandpilebottomisreducedinsoftsoilsubgradereinforcedbychemicalchurningpiles;butthemethodhaslittleeffectontheverticalsettlementofsoilunderpilebottom.Largerpilediameterandhighermaterialstrengthcanimprovethebearingcapacityofcompositesubgrade.ButdifferentlayoutsofchemicalchurningpilesandwhethertosetupGoodmaninterfaceelementhavelittleinfluenceonbearingcapacityofcompositesubgrade.Keywords:chemicalchurningpile;compositesubgrade;bearingcharacteristics;three-dimensionalnumericalanalysis1引言高压旋喷注浆法是将带有特殊喷嘴的注浆管置于土层预定深度,以高压喷射流将固化浆液与土体混合、凝固硬化加固地基的方法[1]。若在喷射的同时,喷嘴以一定的速度旋转、提升,则形成喷浆液与土混合的圆柱形桩体,通常称为旋喷桩[2]。高压旋喷桩地基加固技术在20世纪70年代初发展起来,之后在国内外发展十分迅速。目前,对旋喷桩的研究多数集中在工法的改进上[3-4],或者针对单桩的承载特性进行研究,而对旋喷群桩的承载特性则报道很少。本文根据工程实际情况,采用基于MIDAS-GTS的三维有限元分析技术对旋喷群桩复合地基承载特性进行了研究。2有限元计算模型为了便于分析在旋喷群桩的布置方式、桩弹性模量、桩径及桩距等设计参数及桩-土接触面参数变化下复合地基的承载特性,本文建立了用于对比的基准有限元分析模型,并通过改变基准模型中的对第3期安关峰等:旋喷群桩复合地基承载特性的数值分析应参数进行计算比较得出结论。基准模型中,旋喷桩桩长为10m,桩径为500mm,桩距为1m。土层共2层,其中上层的土层1厚度为6m。复合地基的上部荷载采用均布荷载,数值为90kPa。土层及旋喷桩桩体均采用M-C本构材料模型。模型四周及底部均为对应法向方向的平移约束。基准模型的总体单元数量为98100个,节点数量为51956个,所有单元均为六节点五面体实体单元。基准的有限元整体模型不考虑在桩-土之间设置Goodman接触面单元[5-6](在3.6节中专门阐述了接触面单元设置对计算结果的影响。整体及旋喷桩模型如图1所示,旋喷桩及土层的相关参数如表1所示。(a整体模型(b旋喷桩网格模型图1整体有限元分析模型及旋喷桩网格模型Fig.1Finiteelementmodelsoffoundationandchemicalchurningpiles表1旋喷桩桩体及土层材料参数表Table1Physicalparametersofsoilandchemicalchurningpile层号土类名称弹性模量/MPa泊松比重度/(kN/m3黏聚力/kPa内摩擦角/(°1旋喷桩104000.2021.5900382土层1150.3518.020102土层2400.3119.040203旋喷桩复合地基参数对地基承载特性的影响3.1旋喷桩布置方式对复合地基承载特性的影响图2为旋喷桩矩形布置和梅花形布置示意图。通过分析可知,在本文设定的参数条件下,当旋喷桩布置方式为矩形布置时(见图3,复合地基的最大竖向沉降为12.5mm,旋喷桩桩体的最大竖向应力为634.9kPa;当旋喷桩布置方式为梅花形布置时(见图4,复合地基的最大竖向沉降为12.4mm,旋喷桩桩体的最大竖向应力为637.8kPa(位置均在旋喷桩桩底上方1.5m的桩身位置。从数值上可以看出,在不同的布置方式下,竖向沉降量及旋喷桩竖向位移的差值与绝对数值的比值均在0.5%以内。由此可以得出结论,在旋喷桩置换率一定的情况下,复合地基采用这两种不同的旋喷桩布置方式时,其对复合地基的变形及受力特征影响很小,可以忽略不计。(a矩形布置(b梅花形布置图2旋喷桩布置方式Fig.2Distributionsofchemicalchurningpiles(a地基沉降(b旋喷桩竖向应力图3矩形布置时地基的沉降及旋喷桩的竖向应力分布Fig.3Settlementofsubgradeandverticalstressofchemicalchurningpilewithrectangledistribution(a地基沉降(b旋喷桩竖向应力图4梅花形布置时地基沉降及旋喷桩的竖向应力分布Fig.4Settlementofsubgradeandverticalstressofchemicalchurningpilewithquincunxdistribution3.2旋喷桩桩长对复合地基承载特性的影响本文基准模型中旋喷桩的桩长为10m。为了研沉降/mm-0.8-1.6-2.3-3.1-3.9-4.7-5.4-6.2-7.0-7.8-8.6-9.4-10.1-10.9-11.7-12.5沉降/mm-0.8-1.6-2.3-3.1-3.9-4.6-5.4-6.2-7.0-7.8-8.5-9.3-10.1-10.9-11.6-12.4-174.6-203.5-232.5-261.4-290.4-319.4-348.3-377.3-406.2-435.2-464.1-493.1-522.0-551.0-579.9-608.9-637.8应力/kPa907岩土力学2021年究旋喷桩的桩长对复合地基承载特性的影响,分别建立了桩长为5~19m(以1m为增量的数值分析模型,并将计算结果与基准模型结果进行比较(见图5。由图5可知,在桩长不变的情况下,桩间土在地表的竖向沉降值最大;从地表至下方1m深度内的土体竖向沉降有一定量的减少;从地表下方1m至桩底深度之间的土体沉降值基本不变;而从桩底深度往下竖向沉降值则呈线性递减。在地面荷载一定的情况下,随着旋喷桩桩长的增加,桩间土体地表处的竖向沉降呈线性减少趋势。从图5还可看出,在旋喷桩桩长改变的情况下,桩底下方土体的竖向沉降曲线是基本重合的。这也说明旋喷桩加固复合地基主要是减小了地表至桩底范围内土体的竖向沉降值,而对下方的土体沉降基本无影响。图5不同桩长时桩间土竖向沉降随埋深的变化曲线Fig.5Relationshipsbetweendepthandverticalsettlementofsoilbetweenpileswithdifferentpilelengths由图6可知,在桩长一定的情况下,旋喷桩桩体的竖向应力值在从地表至下方2m深度范围内的增速较大(该段范围内桩间土体竖向沉降比桩体的大,对旋喷桩产生了向下的摩擦力;之后从地表下方2m深度起至旋喷桩桩底上方1m范围内竖向应力值继续增加,但增速减小(该段范围内土体与桩体的竖向位移逐渐趋于一致,二者共同变形。从图6可看出,在旋喷桩桩长改变的情况下,所有旋喷桩从地表至桩底上方1m范围内的竖向应力增加曲线基本重合于同一条曲线。这说明桩长的不同并未改变旋喷桩桩身竖向应力随深度的分布趋势。除了临近桩底的部分以外,不同桩长的旋喷桩桩身竖向应力随深度的增加曲线基本重合于同一条曲线。图6不同桩长时桩身竖向应力随桩身深度的变化曲线Fig.6Relationshipsbetweenverticalstressofpileanddepthwithdifferentpilelengths3.3旋喷桩弹性模量对复合地基承载特性的影响本文中基准模型的弹性模量E=10400MPa。为了研究旋喷桩弹性模量对复合地基承载特性的影响,建立了弹性模量为0.25E、0.50E、0.75E及1.25E、1.50E的数值模型,并将计算结果与基准模型结果进行比较。由图7、8可知,在地表均布荷载的作用下,旋喷桩桩身竖向应力在地表附近(本文中为地表至地表下方1.5m深度范围内迅速增加,到达一定深度后增速减小,并在临近桩底深度之前竖向应力由增变减。相应地,桩间土体竖向应力在地表至下方1.5m范围内呈减小趋势,再往下则是随着深度的增加而增加,其中在旋喷桩桩底附近深度的增速较大。由图7、8可知,旋喷桩弹性模量的变化对桩身竖向应力及桩间土竖向应力的分布趋势影响均不大。图9为旋喷桩弹性模量分别为0.25E、1.00E及1.50E时的复合地基整体竖向应力分布情况。为了研究旋喷桩弹性模量对地基承载特性的影响,选取旋喷桩竖向应力最大位置处(本文为地表7.5m处桩身竖向应力进行研究。图10为该深度处旋喷桩桩身竖向应力值随弹性模量的变化曲线。由图可知,随着旋喷桩弹性模量的增加,同一位置处旋喷桩的竖向应力也在增加,但增速呈减小趋势。这说明提高旋喷桩材料的弹性模量会提高复合地基的承载能力。但旋喷桩材料强度与地基承载能力并不是呈线性关系,当旋喷桩材料弹性模量达到一定值后,继续增加对提高复合地基承载能力的贡献不大。5m6m7m8m9m10m11m12m13m14m15m16m17m18m19m5m6m7m8m9m10m11m12m13m14m15m16m17m18m19m908第3期安关峰等:旋喷群桩复合地基承载特性的数值分析图7不同弹性模量时桩身竖向应力随深度的变化曲线Fig.7Relationshipsbetweenverticalstressofpileanddepthwithdifferentmoduliofelasticity图8不同弹性模量时桩间土体竖向应力随深度的变化曲线Fig.8Relationshipsbetweendepthandverticalstressofsoilbetweenpileswithdifferentmoduliofelasticity(a0.25E(b1.00E(c1.50E图9不同弹性模量时复合地基的竖向应力分布Fig.9Verticalstressesofcompositesubgradewithdifferentmoduliofelasticityofpile图10埋深7.5m处桩身竖向应力值随弹性模量的变化曲线Fig.10Relationshipsbetweenverticalstressandmodulusofelasticityofchemicalchurningpilewithdepthof7.5m3.4旋喷桩桩径对复合地基承载特性的影响本文中基准模型的桩径为0.5m。为了研究旋喷桩桩径对复合地基承载特性的影响,分别建立了桩径为0.4、0.6、0.8m的数值分析模型(桩间距均为1m,并将计算结果与基准模型结果进行比较。由图11可知,在上部荷载一定的情况下,随着桩径的增加,桩身内部的竖向应力随之降低。这是由于在桩径较小时,由于复合地基置换率低,桩-土之间弹性模量的差异导致了旋喷桩承受了上部荷载的绝大部分,因此,桩身应力比较高。由于桩径的增加,旋喷桩复合地基的整体置换率提高,更多比率的旋喷桩桩体参与承担上部荷载,因此,相对而言,桩身应力就降低。进一步分析可知,在旋喷桩桩距不变的情况下,随着桩径的增加,旋喷桩复合地基的承载能力也随之提高。由图12可知,地基在地表处的沉降值比较大。对于某一桩径时的复合地基而言,从地表开始至下方较小深度范围内(本文中该范围为0~-2m,地基沉降值迅速减小。随着桩径的增加,从地表开始至下方较小深度范围内的地基沉降值会逐渐减小;而再往下(本算例中为地表下方-2m深度以下的地基沉降值则基本不随桩径增加而变化。可见,桩径的改变主要影响地表至下方较小深度范围内的地基竖向沉降。应力/kPa-37.2-72.6-108.1-143.6-179.0-214.5-250.0-285.4-320.9-356.4-391.8-427.3-462.8-498.3-533.7-539.2-604.7应力/kPa634.6应力/kPa-34.4-72.2-109.9-147.6-185.4-223.1-260.8-298.6-336.3-374.1-411.8-449.5-487.3-525.0-562.7-600.5-638.2桩身标高/m桩身竖向应力/kPa0-1-2-3-4-5-6-7-8-9-10100200300400土层标高/m0.25E0.50E0.75E1.00E1.25E1.50E0-2-4-6-8-10-12-14-16土体竖向应力/kPa909岩土力学2021年图11不同桩径时桩身竖向应力值随桩身深度的变化曲线Fig.11Relationshipsbetweenverticalstressanddepthwithdifferentpilediameters图12不同桩径时桩间土竖向沉降值随深度的变化曲线Fig.12Relationshipsbetweendepthandverticalsettlementofsoilbetweenpileswithdifferentpilediameters为了研究随着旋喷桩桩径增加对复合地基承载特性的影响,选取桩间土的地表沉降进行研究。图13为桩间土的地表沉降值随桩径的变化曲线。由图可知,随着旋喷桩桩径的增加,桩间土的地表沉降随之减小,但减速呈降低的趋势。这说明增加旋喷桩桩径会减少复合地基土地表的竖向沉降。但旋喷桩桩径的增加与地基抵抗竖向沉降的能力并不是呈线性关系,当旋喷桩桩径达到一定值后,继续增大桩径对提高地基抵抗竖向沉降能力的贡献不大。图13桩间土竖向沉降值随桩径的变化曲线Fig.13Relationshipbetweenverticalsettlementanddiameterofchemicalchurningpile3.5旋喷桩桩距对复合地基承载特性的影响本文中基准模型的桩距为1m。为了研究旋喷桩桩距对复合地基承载特性的影响,分别建立了桩径为0.5m,桩距为0.6、0.8、1.5m的数值分析模型,并将计算结果与基准模型结果进行比较。由图14可知,在上部荷载一定的情况下,随着桩距的增加,桩身内部的竖向应力随之增加。这是由于在桩距较小时,由于复合地基置换率比较高,旋喷桩桩体的竖向应力比较低。随着桩距的增加,旋喷桩复合地基的整体置换率降低,旋喷桩承担了更多的上部荷载,因此,相对而言,桩身应力提高。当桩距增加到一定距离后,桩体本身的竖向应力会进一步提高,直至达到抗压强度而发生破坏。图14不同桩距时桩身竖向应力值随深度的变化曲线Fig.14Relationshipsbetweenverticalstressanddepthwithdifferentpiledistances43211-0.4m2-0.5m3-0.6m4-0.8m1-0.4m2-0.5m3-0.6m4-0.8m123443211-1.5m2-1.0m3-0.8m4-0.6m910第3期安关峰等：旋喷群桩复合地基承载特性的数值分析911由图15可知，在桩距较小情况下，地基在地表处的沉降值比较小。而随着桩距的增加，地基在地表处的沉降值逐渐增大。特别是当桩距增加到一定程度时（如本文算例中桩距为1.5m时），不仅在地表处的桩间土竖向沉降较大，在地表下方一定深度内的竖向位移值较其他桩距（桩距为0.6、0.8、1m时）也更大。这说明当桩距增大至一定数值时，旋喷桩同桩间土之间失去了整体变形协调的能力，复合地基的承载能力也大为降低。3.6桩-土界面单元对复合地基承载特性的影响为了研究在旋喷桩桩体与周边土体之间设置Goodman接触面单元对复合地基承载特性的影响，本文建立了考虑桩-土接触的复合地基数值分析模型。其中接触面参数的取值按照旋喷桩桩体以及周边土体的特性并结合以往同类文献的经验进行选定，接触面单元的法向刚度取值为5×106kN/m3，切向刚度的取值为5×105kN/m3。图17为整体有限元模型中的接触面单元网格。43211－1.5m2－1.0m3－0.8m4－0.6m图17接触面单元网格Fig.17Interfaceelementmeshes图18、19分别为考虑桩-土接触以及未考虑接触时复合地基的沉降及竖向应力分布。由图可知，在考虑接触面单元时，复合地基的最大竖向沉降为图15不同桩距时桩身竖向沉降值随深度的变化曲线Fig.15Relationshipsbetweenverticalsettlementanddepthwithdifferenepiledistances12.6mm，旋喷桩桩体的最大竖向应力为632.6kPa；未考虑接触面单元时，复合地基的最大竖向沉降为12.5mm，旋喷桩桩体的最大竖向应力为634.6kPa。通过计算发现，当桩-土之间设置与未设置Goodman接触面单元时相比地基的受力及变形情况均差别不大。但设置了接触面单元的情况下，旋喷桩桩周土体的竖向应力受桩体影响的区域相较未设接触面时会相对小一些。沉降/mm0-0.8-1.6-2.4-3.2-3.9-4.7-5.5-6.3-7.1-7.9-8.7-9.4-10.2-11.0-11.8-12.6应力/kPa-33.3-70.8-108.2-145.7-183.1-220.6-258.0-295.5-332.9-370.4-407.8-445.3-482.8-520.2-557.7-595.1-632.6为了研究随着旋喷桩桩距增加对复合地基承载特性的影响，选取桩间土的地表沉降进行研究。图16为桩间土的地表沉降值随桩距的变化曲线。由图可知，随着旋喷桩桩距的增加，桩间土的地表沉降随之增加，且增速呈加快的趋势。这说明增加旋喷桩桩距会加速增加复合地基的竖向沉降，旋喷桩复合地基承载能力加速降低。(a地基沉降(b地基竖向应力图16桩间土竖向沉降值随桩距的变化曲线Fig.16Relationshipbetweenverticalsettlementanddistancebetweenpiles图18考虑接触时地基的沉降及竖向应力分布（剖切图）Fig.18Verticalsettlementandstressofsubgradeconsideringinterfaceelement912沉降/mm0-0.8-1.6-2.3-3.1-3.9-4.7-5.4-6.2-7.0-7.8-8.6-9.4-10.1-10.9-11.7-12.5岩应力/kPa-34.7-72.2-109.7-147.2-184.7-222.2-259.6-297.1-334.6-372.1-409.6-447.1-484.6-522.1-559.6-597.1-634.6土力学2021年参考文献[1]姚贤华,裴松伟,赵顺波.高压旋喷桩复合地基承载特性的有限元分析[J].华北水利水电学院学报,2021,30(1:93－95.YAOXian-hua,PEISong-wei,ZHAOShun-bo.Finiteelementanalysisofload-carryingcapacitywithcompositefoundationofhighpressurerotarygroutingpile[J].JournalofNorthChinaInstituteofWater(a地基沉降(b地基竖向应力图19未考虑接触时地基的沉降及竖向应力分布（剖切图）Fig.19VerticalsettlementandstressofsubgradewithoutconsideringinterfaceelementConservancyandHydroelectricPower,2021,30(1:93－95.[2]《地基处理手册》编写委员会.地基处理手册(第二版[M].北京:中国建筑工业出版社,2000.[3]郝峰.高压旋喷桩复合土钉墙Plaxis有限元分析[J].探矿工程(岩土掘进工程,2021,36(9:52－55.HAOFeng.Plaxisfiniteelementanalysisofsupportingstructurewithhigh-pressurejetgroutingpileandcompositesoil-nailingwall[J].ExplorationEngineering(Rock&SoilDrilling&Tunneling,2021,36(9:52－55.[4]朱晞,王根会.铁路桥梁旋喷桩复合地基的三维弹性有限元分析[J].铁道学报,1996,18(6:95－99.ZHUXi,WANGGen-hui.3-DFEManalysisofcompositefoundationstrengthenedbywhirlysprayedcementpileofrailwaybridge[J].JournalofTheChinaRailwaySociety,1996,18(6:95－99.[5]许宏发,吴华杰,郭少平,等.桩土接触面单元参数分析[J].探矿工程,2002,5:10－12.XUHong-fa,WUHua-jie,GUOShao-ping,etal.Studyoftheparametersofpilesoilcontactsurfaceelement[J].ExplorationEngineering,2002,5:10－12.[6]钱晓丽,陶龙光,刘波.竖向载荷作用下单桩接触面性能分析[J].辽宁工程技术大学学报,2007,26(1:59－61.QIANXiao-li,TAOLong-guang,LIUBo.Performanceanalysisofsinglepile-soilinterfaceunderverticalload[J].JournalofLiaoningTechnicalUniversity,2007,26(1:59－61.4结论（1）当旋喷桩复合地基采用不同的旋喷桩布置方式时，其对复合地基的变形及受力特征影响很小，可以忽略不计。（2）在地面荷载一定的情况下，随着旋喷桩桩长的增加，桩间土体地表处的竖向沉降呈线性减小趋势。旋喷桩加固复合地基主要是减小了地表至桩底深度范围内土体的竖向沉降值，而对桩底下方的土体沉降基本无影响。除了桩底区域以外，不同桩长旋喷桩的桩身竖向应力随深度的增加曲线基本重合于同一条曲线。（3）提高旋喷桩材料的弹性模量会提高复合地基的承载能力，当旋喷桩材料弹性模量达到一定值后，继续增加对提高复合地基承载能力的贡献不大。（4）随着旋喷桩桩径的增加，桩间土的地表沉降随之减少，但减速呈降低的趋势。当旋喷桩桩径达到一定值后，继续增大桩径对提高地基抵抗竖向沉降能力的贡献不大。（5）随着桩距的增加，旋喷桩桩身应力提高。当桩距增加到一定距离后，桩体本身的竖向应力会进一步提高直至达到抗压强度而发生破坏。增加旋喷桩桩距会加速增大复合地基的竖向沉降，旋喷桩复合地基承载能力加速降低。（6）桩-土之间设置与未设置Goodman接触面单元时相比，地基的受力及变形情况均差别不大，但桩周土体竖向应力受影响的区域相较未设接触面时会相对小一些。矩阵的表示理论及其在数值计算中的应用【摘要】：矩阵的理论和方法不仅是各数学学科的基本工具，而且在理论物理学、经济学、统计学、最优化、信息处理、自动控制、工程技术和运筹学等应用学科的理论研究和数值计算中都有着广泛的应用。近年来，随着近代量子力学的不断发展，力学工作者遇到并提出了一系列有关矩阵的理论和计算方面的疑难问题，这些问题制约着量子力学的发展，急需数学工作者给以解答。在本文中，我们通过引入复矩阵的实表示、四元数矩阵的复表示、友向量和伴向量的方法，研究并解决了量子力学等学科中的有关矩阵理论与计算中的下列三类系列疑难问题：1．矩阵的合相似问题两个复矩阵A，B称为是合相似的是指存在复可逆矩阵S满足S~(-1)A(?)=B。我们通过复矩阵的实表示、友向量和伴向量方法，研究并解决了合相似意义下矩阵的若当标准形、合相似意义下矩阵的三角化和矩阵的广义对角化的问题。不但从理论上给出复矩阵的一个新的若当标准形，而且还给出了相应复矩阵若当标准形的计算方法。进一步地，我们不但给出求一个复矩阵A的合若当标准形J的简单方法，而且还给出一种求相应合相似可逆矩阵S(满足S~(-1)A(?)=J)的算法。2．矩阵方程的解问题矩阵方程解的问题是矩阵理论中的一类重要的问题。如何给出某个矩