土木工程中无应力时长度计算程序,可按内容适当作修改再应用于不同的地方

syms P1 P2
%输入参数
format short
syms P1 P2                                %P1-水平力 P2-竖向力
A=[5.0128 6.43384 12.3538]*10^(-3);                          %单位：m2
E=1.95*10^11;                            %单位：N/m2
%Lu= 158.0245;                                %单位：m
w=[407 518 987];
Lx0=[36.596 83.751 127.506]
Ly0=[58.994 66.153 73.697]
FLu0(1,:)=[-15.515 -58.684 -111.569 -164.344 -533.437 -797.016 -1587.734]*1000;  %单位：N
FLu0(2,:)=[48.403 169.178 405.304 640.909 1190.370 1582.780 3937.018]*(-1000);  %单位：N
FLu0(3,:)=[195.294 284.834 460.436 1328.802 2195.121 3061.148 787254.400]*(-1000);  %单位：N
for n=2:2
    Lc=sqrt(Lx0(n)^2+Ly0(n)^2)
Lu1=Lc 
Lu2=1.001*Lc
    for m=4:4
Lu0=[Lu1 Lu2]'
P3=-P1;
P4=w(n)*Lu0-P2;
Ti=sqrt(P1^2+P2^2);
Tj=sqrt(P3^2+P4.^2);
Ly=1/(2*E*A(n)*w(n))*(Tj.^2-Ti^2)+(Tj-Ti)/w(n);
Lx=-P1*(Lu0/E/A(n)+1/w(n)*log((P4+Tj)/(Ti-P2)));
for i=1:2
[FLu(i) q2]=solve(Lx(i)-Lx0(n),Ly(i)-Ly0(n))
end
Lu(n,m)=Lu0(2)+(FLu0(n,m)-FLu(2))*(Lu0(2)-Lu0(1))/(FLu(2)-FLu(1))
Lu0=[Lu0(2) Lu(n,m)]
FLu=[FLu(2) FLu0(n,m)]
for j=1:6
P3=-P1;
P4=w(n)*Lu0-P2;
Ti=sqrt(P1^2+P2^2);
Tj=sqrt(P3^2+P4.^2);
Ly=1/(2*E*A(n)*w(n))*(Tj.^2-Ti^2)+(Tj-Ti)/w(n);
Lx=-P1*(Lu0/E/A(n)+1/w(n)*log((P4+Tj)/(Ti-P2)));
for i=1:2
[FLu(i) q2]=solve(Lx(i)-Lx0(n),Ly(i)-Ly0(n))
end
Lu(n,m)=Lu0(2)+(FLu0(n,m)-FLu(2))*(Lu0(2)-Lu0(1))/(FLu(2)-FLu(1))
Lu0=[Lu0(2) Lu(n,m)]
FLu=[FLu(2) FLu0(n,m)]
end
Lu
end
end