4楼: Originally posted by 风起天阑1 at 2018-09-18 14:10:21
下面图片就是最后数据图的结果，楼主现在只是初学者，只能拾人牙慧，简单的编程，对利用矩阵等方法还不是太过熟悉，还希望大佬不吝赐教，能否帮忙利用矩阵的方法修改一下 。谢谢


...

clear

X = zeros(121,1); % 坐标节点数组
H = zeros(121,1); % 膜厚节点数组
P = zeros(121,1); % 压力节点数组
H1=1.0; % 量纲化最大膜厚
H2=0.5; % 量纲化最小膜厚
U=1;    % 量纲化速度
EDA=0.02;% 润滑油粘度
AL=0.01; % 滑块长度
KG=1;% 等于1时为直线滑块，2时为曲线滑块
N=121;% 节点取值
if KG==2
    X1=-1.0;
    X2=1.0;
else
    X1=0.0;
    X2=1.0;
end
DX=1/(N-1);
if KG == 1 % 直线滑块膜厚计算
    X = X1 : ( X2 - X1 ) / ( N - 1 ) : X2;
    H = H1 / H2 - ( H1 / H2 - 1.0 ) * X;
else     % 曲线滑块膜厚计算
    X = X1 : ( X2 - X1 ) /  ( N - 1 ) : X2;
    H = 1.0 + ( H1 / H2 - 1.0 ) * X .^ 2;
end
P( 1 ) = 0.0; % 压力起点初始化
P( 2 : N -1 ) = 0.5; % 压力各点赋初值
P( N ) = 0.0; % 压力终点初始化
IK=0;
ERO=1;
while ERO>1e-7
    C1=0.0;
    ALOAD=0.0;
    for i=2:N-1
    A1=(0.5*(H(i+1)+H(i)))^3;
    A2=(0.5*(H(i)+H(i-1)))^3;
    PD=P(i);
    P(i)=(-0.5*DX*(H(i+1)-H(i-1))+A1*P(i+1)+A2*P(i-1))/(A1+A2);
    P2(i)=0.3*PD+0.7*P(i);
    if P2(i)<0
       P2(i)=0;
    end
        C1=C1+abs(P2(i)-PD);
        ALOAD=ALOAD+P2(i);
    end
   ERO=C1/ALOAD; % 相对误差，用于判断收敛情况
   IK=IK+1;
end

X0 = sum( P .* X' ) * AL;
ALOAD=ALOAD*DX*AL*6.0*U*EDA*AL/(H2^2);% 承载能力计算
P0 = zeros(121,1);
if KG==1
     for i=1:N
       P0( i )=-(-1.0/(H(i)*H2)+H1*H2/(H1+H2)/(H2*H(i))^2+1.0/(H1+H2))/(H1/H2-1.0)*H2;
     end
end
figure
subplot( 2, 2, 1 )
plot( X, P0, '--b' )
subplot( 2, 2, 3 )
plot( X, H, 'r' )
subplot( 2, 2, [ 2, 4 ] )
plot( X, P0, '--b' )
hold on
plot( X, H, 'r' )

5楼: Originally posted by xiegangmai at 2018-09-18 15:27:18
改了改部分代码，供参考：

clear

X = zeros(121,1); % 坐标节点数组
H = zeros(121,1); % 膜厚节点数组
P = zeros(121,1); % 压力节点数组
H1=1.0; % 量纲化最大膜厚
H2=0.5; % 量纲化最小膜厚
U=1;    ...

4楼: Originally posted by 风起天阑1 at 2018-09-18 14:10:21
下面图片就是最后数据图的结果，楼主现在只是初学者，只能拾人牙慧，简单的编程，对利用矩阵等方法还不是太过熟悉，还希望大佬不吝赐教，能否帮忙利用矩阵的方法修改一下 。谢谢


...