源程序.txt

本软件研究了地基的计算方法

%......没用迭代计算扰度的程序..........%
%......D为变值,g为改进过..........%
function Calculate_g(row,column,long,wide,E,h,p,coefficient,force,a,EI1,EI2,EI3,EI4)
% 输入参数，row,column为横纵方向上的节点数，长，宽，弹性模量，板厚，泊松比，地基系数,集中力
h1=long/(row-1);%X方向上的步长
h2=wide/(column-1);%Y方向上的步长
t1=E*h^3/(12*(1-p^2));%暂存值
D=zeros(row,column);%赋值
for i=1:row
    for j=1:column
        D(i,j)=t1*(1-4*(1-a)*(((i-1)*h1/long-1/2)^2+((j-1)*h2/wide-1/2)^2));
    end
end
DL=zeros(row-1,column-1);%赋值
ul=zeros(row-1,column-1);%赋值
for i=1:(row-1)
    for j=1:(column-1)
        DL(i,j)=t1*(1-4*(1-a)*(((i-0.5)*h1/long-1/2)^2+((j-0.5)*h2/wide-1/2)^2));
        ul(i,j)=p;
    end
end
u=zeros(row,column);
C=zeros(row,column);
for i=1:row
    for j=1:column
        C(i,j)=coefficient;%地基系数
        u(i,j)=p;
    end
end
t2=force/(h1*h2);%计算作用力
count=row*column;%count为方程的系数矩阵的维数，行*列%
%...................调用扰度函数.................................. %
DiJi(row,column,h1,h2,D,DL,u,ul,t2,C,EI1,EI2,EI3,EI4,count)
%梁1,2,3,4的取值，初值设为0，常数 %
%........................................................
%计算扰度的函数,A为系数矩阵,常数项b,X为扰度,row行节点数,column为列节点数,h1,h2为步长,u为常数
%D,u，c分别为各点的弯曲刚度值与...，需要输入size(D)=row*column维矩阵，u，c维数也一样
%DL为节点之中间的弯曲刚度值，输入size(DL)=row*column，ul也一样
%q为输入值，size(q)=row*column,1%
%                row=N1+1，column=N2+1          %
function DiJi(row,column,h1,h2,D,DL,u,ul,t2,C,EI1,EI2,EI3,EI4,count)
%D代表刚度系数，DL代表/D的刚度系数%
%    计算方程的系数矩阵A，标号从(1,1)算起  %
bc0=4;%A对角线上的步长%
bc1=1;%A对角线上一行的步长%
bc2=1;%A对角线上两行的步长%
bc3=4;%S对角线的步长%
bc4=1;%S对角线上一行的步长%
bc5=1;%S对角线下一行的步长%
bc6=4;%T对角线的步长%
h12=h1^2*h2^2;
k1=(h2/h1)^2;
k2=1/k1;
N1=row-1;%N1为地基最右端的坐标,N1>=7，行坐标
N2=column-1;%N2为地基最上端的坐标,N2>=7，列坐标
q=zeros(count,1);%常数项                                                                                                      
b=zeros(count,1);%常数项
%.........地基为N1*N2 ....需改动A,B,C,D点的位置...
q(column*ceil(N1/4)+ceil(column/4),1)=t2;%A点作用力,X1，X2坐标为长(N1),高(N2)的的四分子一，
q(column*ceil(N1*3/4)+ceil(column/4),1)=t2;%B点作用力，X1坐标为长(N1),高(N2)的四分子三
q(column*ceil(N1/4)+ceil(column*3/4),1)=t2;%C点作用力，X2坐标为ceil(N1/4)=5
q(column*ceil(N1*3/4)+ceil(column*3/4),1)=t2;%D点作用力，X1，X2坐标为长(N1),高(N2)的的四分子三
%.............对q,b稀疏化，节省内存.........................
q=sparse(q);
b=sparse(b);
%.......................................................
A=zeros(count);%生成n*n维的零矩阵，以便对其赋值%
%.............对A稀疏化，节省内存........................
A=sparse(A);
%......................................................
%原文中的DL(0.5,0.5)为DL(1,1)勿需改动,D(1,0)为D(2,1),u(0,1)为u(1,2),ul与原文相同，C(0,0)为C(1,1)，循环时DL，ul跟D，C都得变化,但是x坐标不变，y坐标变时不加1，D，C的x，y坐标都在原来基础上加1%
%                 输入A中的数据元素：对角线上一行的元素         %
%             下面为地基中的第一列               %
  A(1,2)=-2*DL(1,1)*(1-ul(1,1))-k2*(D(1,2)*(1-u(1,2)*u(1,2))+2*EI1/h1)-5*EI1*k2/h1;
%a'(0,1)=-2*DL(1,1)*(1-ul(1,1))-k2*(D(0,1)*(1-u(0,1)*u(0,1))+2*EI1/h1)-5*EI1*k2/h1%
%输入纸1，/a01,W01的系数%
%点(0,1)的D(0,1)相当于矩阵中的D(1,2)%
  A(2,3)=-2*DL(1,2)*(1-ul(1,2))-k2*(D(1,2)*(1-u(1,2)*u(1,2))+2*EI1/h1)-k2*(D(1,3)*(1-u(1,3)*u(1,3))+2*EI1/h1)-10*EI1*k2/h1;
%a'(0,2)=-2*DL(1,2)*(1-ul(1,2))-k2*(D(0,1)*(1-u(0,1)*u(0,1))+2*EI1/h1)-k2*(D(0,2)*(1-u(0,2)*u(0,2))+2*EI1/h1)-10*EI1*k2/h1%
%输入纸2，/a02,W02的系数%
  A(3,4)=-2*DL(1,3)*(1-ul(1,3))-k2*(D(1,3)*(1-u(1,3)*u(1,3))+2*EI1/h1)-k2*(D(1,4)*(1-u(1,4)*u(1,4))+2*EI1/h1)-2*EI1*k2/h1;
%a'(0,3)=-2*DL(1,3)*(1-ul(1,3))-k2*(D(0,2)*(1-u(0,2)*u(0,2))+2*EI1/h1)-k2*(D(0,3)*(1-u(0,3)*u(0,3))+2*EI1/h1)-2*EI1*k2/h1%
%输入纸3，/a03,W03的系数%
  A(4,5)=-2*DL(1,4)*(1-ul(1,4))-k2*(D(1,4)*(1-u(1,4)*u(1,4))+2*EI1/h1)-k2*(D(1,5)*(1-u(1,5)*u(1,5))+2*EI1/h1);
%a'(0,4)=-2*DL(1,4)*(1-ul(1,4))-k2*(D(0,3)*(1-u(0,3)*u(0,3))+2*EI1/h1)-k2*(D(0,4)*(1-u(0,4)*u(0,4))+2*EI1/h1)%
%输入纸4，/a04,W04的系数%
%N2为列节点末点的坐标，通式%
for j=5:N2-3  
    A(j,j+1)=-2*DL(1,j)*(1-ul(1,j))-k2*(D(1,j)*(1-u(1,j)*u(1,j))+2*EI1/h1)-k2*(D(1,j+1)*(1-u(1,j+1)*u(1,j+1))+2*EI1/h1);
    %a'(0,j)=-2*DL(1,j)*(1-ul(1,j))-k2*(D(0,j-1)*(1-u(0,j-1)*u(0,j-1))+2*EI1/h1)-k2*(D(0,j)*(1-u(0,j)*u(0,j))+2*EI1/h1)%
    %输入纸5，/a0k,W0(j+1)的系数,其中j=5,6,...,N2-3%
end
A(N2-2,N2-1)=-2*DL(1,N2-2)*(1-ul(1,N2-2))-k2*(D(1,N2-2)*(1-u(1,N2-2)*u(1,N2-2))+2*EI1/h1)-k2*(D(1,N2-1)*(1-u(1,N2-1)*u(1,N2-1))+2*EI1/h1)-2*EI1*k2/h1;
% a'(0,N2-2)=-2*DL(1,N2-2)*(1-ul(1,N2-2))-k2*(D(0,N2-3)*(1-u(0,N2-3)*u(0,N2-3))+2*EI1/h1)-k2*(D(0,N2-2)*(1-u(0,N2-2)*u(0,N2-2))+2*EI1/h1)-2*EI1*k2/h1%
%输入纸6，/a0(N2-2),W0(N2-2)的系数%
 A(N2-1,N2)=-2*DL(1,N2-1)*(1-ul(1,N2-1))-k2*(D(1,N2-1)*(1-u(1,N2-1)*u(1,N2-1))+2*EI1/h1)-k2*(D(1,N2)*(1-u(1,N2)*u(1,N2))+2*EI1/h1)-10*EI1*k2/h1;
%a'(0,N2-1)=-2*DL(1,N2-1)*(1-ul(1,N2-1))-k2*(D(0,N2-2)*(1-u(0,N2-2)*u(0,N2-2))+2*EI1/h1)-k2*(D(0,N2-1)*(1-u(0,N2-1)*u(0,N2-1))+2*EI1/h1)-10*EI1*k2/h1%
%输入纸7，/a0(N2-1),W0(N2-1)的系数%
A(N2,column)=-2*DL(1,N2)*(1-ul(1,N2))-k2*(D(1,N2)*(1-u(1,N2)*u(1,N2))+2*EI1/h1)-5*EI1*k2/h1;
% a'(0,N2)=-2*DL(1,N2)*(1-ul(1,N2))-k2*(D(0,N2-1)*(1-u(0,N2-1)*u(0,N2-1))+2*EI1/h1)-5*EI1*k2/h1%
%输入纸8，/a0N2,W0N2的系数%

%             下面为地基中的第二列到N1列 i=1,2...N1-1               %
for i=N2+2:column:((N1-1)*column+1)          %从N2+2处开始,步长为N2+1%
    A(i,i+1)=-2*D(i-bc1*N2,2)*(k2+u(i-bc1*N2,2))-2*DL(i-bc1*N2-1,1)*(1-ul(i-bc1*N2-1,1))-2*DL(i-bc1*N2,1)*(1-ul(i-bc1*N2,1));
    %a'(i,1)=-2*D(i,1)*(k2+u(i,1))              -2*DL(i,1)*(1-ul(i,1))-2*DL(i+1,1)*(1-ul(i+1,1))%
    %输入论文55，/ai.1%
    for j=i+1:i+N2-2
        A(j,j+1)=-2*D(i-bc1*N2,j-i+1)*(k2+u(i-bc1*N2,j-i+1))-2*D(i-bc1*N2,j-i+2)*(k2+u(i-bc1*N2,j-i+2))-2*DL(i-bc1*N2-1,j-i+1)*(1-ul(i-bc1*N2-1,j-i+1))-2*DL(i-bc1*N2,j-i+1)*(1-ul(i-bc1*N2,j-i+1));
      %  a'(i,j)=-2*D(i,j-1)*(k2+u(i,j-1))                  -2*D(i,j)*(k2+u(i,j))                      -2*DL(i,j)*(1-ul(i,j))                          -2*DL(i+1,j)*(1-ul(i+1,j))%
        %输入论文55,/aij，其中j=2,...,N2-1%
    end
    A(i+N2-1,i+N2)=-2*D(i-bc1*N2,N2)*(k2+u(i-bc1*N2,N2))-2*DL(i-bc1*N2-1,N2)*(1-ul(i-bc1*N2-1,N2))-2*DL(i-bc1*N2,N2)*(1-ul(i-bc1*N2,N2));
   %      a'(i,N2)=-2*D(i,N2-1)*(k2+u(i,N2-1))        -2*DL(i,N2)*(1-ul(i,N2))                -2*DL(i+1,N2)*(1-ul(i+1,N2))%
    %输入论文55,/aiN2%
    bc1=bc1+1;
end

%             下面为地基中的第N1+1列               %
  A(N1*column+1,N1*column+2)=-2*DL(N1,1)*(1-ul(N1,1))-k2*(D(row,2)*(1-u(row,2)*u(row,2))+2*EI3/h1)-5*EI3*k2/h1;
%                   a'(N1,1)=-2*DL(N1,1)*(1-ul(N1,1))-k2*(D(N1,1)*(1-u(N1,1)*u(N1,1))+2*EI3/h1)-5*EI3*k2/h1%
%输入纸24，/aN1.1,WN1.1的系数%
A(N1*column+2,N1*column+3)=-2*DL(N1,2)*(1-ul(N1,2))-k2*(D(row,2)*(1-u(row,2)*u(row,2))+2*EI3/h1)-k2*(D(row,3)*(1-u(row,3)*u(row,3))+2*EI3/h1)-10*EI3*k2/h1;
%                 a'(N1,2)=-2*DL(N1,2)*(1-ul(N1,2))-k2*(D(N1,1)*(1-u(N1,1)*u(N1,1))+2*EI3/h1)      -k2*(D(N1,2)*(1-u(N1,2)*u(N1,2))+2*EI3/h1)-10*EI3*k2/h1%
%输入纸25，/aN1.2,WN1.2的系数%
A(N1*column+3,N1*column+4)=-2*DL(N1,3)*(1-ul(N1,3))-k2*(D(row,3)*(1-u(row,3)*u(row,3))+2*EI3/h1)-k2*(D(row,4)*(1-u(row,4)*u(row,4))+2*EI3/h1)-2*EI3*k2/h1;
%                 a'(N1,3)=-2*DL(N1,3)*(1-ul(N1,3))-k2*(D(N1,2)*(1-u(N1,2)*u(N1,2))+2*EI3/h1)-k2*(D(N1,3)*(1-u(N1,3)*u(N1,3))+2*EI3/h1)-2*EI3*k2/h1%
%输入纸26，/aN1.3,WN1.3的系数%
A(N1*column+4,N1*column+5)=-2*DL(N1,4)*(1-ul(N1,4))-k2*(D(row,4)*(1-u(row,4)*u(row,4))+2*EI3/h1)-k2*(D(row,5)*(1-u(row,5)*u(row,5))+2*EI3/h1);
%                 a'(N1,4)=-2*DL(N1,4)*(1-ul(N1,4))-k2*(D(N1,3)*(1-u(N1,3)*u(N1,3))+2*EI3/h1)-k2*(D(N1,4)*(1-u(N1,4)*u(N1,4))+2*EI3/h1)%
%输入纸27，/aN1.4,WN1.4的系数%
for j=N1*column+5:row*column-4
    A(j,j+1)=-2*DL(N1,j-N1*column)*(1-ul(N1,j-N1*column))-k2*(D(row,j-N1*column)*(1-u(row,j-N1*column)*u(row,j-N1*column))+2*EI3/h1)-k2*(D(row,j-N1*column+1)*(1-u(row,j-N1*column+1)*u(row,j-N1*column+1))+2*EI3/h1);
    %a'(N1,j)=-2*DL(N1,j)*(1-ul(N1,j))                   -k2*(D(N1,j-1)*(1-u(N1,j-1)*u(N1,j-1))+2*EI3/h1)                              -k2*(D(N1,j)*(1-u(N1,j)*u(N1,j))+2*EI3/h1)%
    %输入纸28,/aN1.k，WN1.(j+1)的系数,其中K=5,6,...,N2-3%
end
A(row*column-3,row*column-2)=-2*DL(N1,N2-2)*(1-ul(N1,N2-2))-k2*(D(row,N2-2)*(1-u(row,N2-2)*u(row,N2-2))+2*EI3/h1)-k2*(D(row,N2-1)*(1-u(row,N2-1)*u(row,N2-1))+2*EI3/h1)-2*EI3*k2/h1;
%                      a'(N1,N2-2)=-2*DL(N1,N2-2)*(1-ul(N1,N2-2))-k2*(D(N1,N2-3)*(1-u(N1,N2-3)*u(N1,N2-3))+2*EI3/h1)      -k2*(D(N1,N2-2)*(1-u(N1,N2-2)*u(N1,N2-2))+2*EI3/h1)-2*EI3*k2/h1%
%输入纸29，/aN1.(N2-2),WN1.(N2-2)的系数%
A(row*column-2,row*column-1)=-2*DL(N1,N2-1)*(1-ul(N1,N2-1))-k2*(D(row,N2-1)*(1-u(row,N2-1)*u(row,N2-1))+2*EI3/h1)-k2*(D(row,N2)*(1-u(row,N2)*u(row,N2))+2*EI3/h1)-10*EI3*k2/h1;
%                      a'(N1,N2-1)=-2*DL(N1,N2-1)*(1-ul(N1,N2-1))-k2*(D(N1,N2-2)*(1-u(N1,N2-2)*u(N1,N2-2))+2*EI3/h1)      -k2*(D(N1,N2-1)*(1-u(N1,N2-1)*u(N1,N2-1))+2*EI3/h1)-10*EI3*k2/h1 %
%输入纸30，/aN1.(N2-1),WN1.(N2-1)的系数%
A(row*column-1,row*column)=-2*DL(N1,N2)*(1-ul(N1,N2))-k2*(D(row,N2)*(1-u(row,N2)*u(row,N2))+2*EI3/h1)-5*EI3*k2/h1;
%                      a'(N1,N2)=-2*DL(N1,N2)*(1-ul(N1,N2))-k2*(D(N1,N2-1)*(1-u(N1,N2-1)*u(N1,N2-1))+2*EI3/h1)-5*EI3*k2/h1 %
%输入纸31，/aN1.N2,WN1.N2的系数%

 %                 输入A中的数据元素：对角线上两行的元素         %
 %             下面为地基中的第一列               %
  A(1,3)=1/2*k2*(D(1,2)*(1-u(1,2)*u(1,2))+2*EI1/h1)+4*EI1*k2/h1;
%a"(0,1)=1/2*k2*(D(0,1)*(1-u(0,1)*u(0,1))+2*EI1/h1)+4*EI1*k2/h1%
%输入纸1，//a0.1,W0.2的系数%
  A(2,4)=1/2*k2*(D(1,3)*(1-u(1,3)*u(1,3))+2*EI1/h1)+5/2*EI1*k2/h1;
%a"(0,2)=1/2*k2*(D(0,2)*(1-u(0,2)*u(0,2))+2*EI1/h1)+5/2*EI1*k2/h1%
%输入纸2，//a0.2,W0.3的系数%
  A(3,5)=1/2*k2*(D(1,4)*(1-u(1,4)*u(1,4))+2*EI1/h1);
%a"(0,3)=1/2*k2*(D(0,3)*(1-u(0,3)*u(0,3))+2*EI1/h1) %