2016年全国数模A题的问题1的MATLAB代码
代码如下：
%2016.9.11 系泊系统设计
% r,H,v,V,beta,rou,L
% r 是钢球质量比（以 1200kg 为基准），即钢球质量=1200kg*r，调整 r 即可调整钢球质量。H 为水深，v 为风速， V为水速，beta 为风速水速之间的夹角，rou 是锚链密度，L 是锚链的长度
%phi0,phi1,phi2,phi3,phi4,h,theta,x
%phi0 是钢桶的倾斜角，phi1 至 phi4 为从下至上各节钢管的倾斜角，h 为吃水深度，theta为锚链末端与锚的连接 点处切线方向与海床夹角，x 为游动的最大范围。
r=1;
H=18;
v=12;
V=0;
beta=0;
rou=7;
L=22.05;
%考虑到锚链可能并不是全部悬空，有可能部分堆放在海床上，因此先求解不堆放情形，至于堆放情形，在下面 if 中考虑。
%fzero函数是用于求解单变量非线性方程的，它需要提供一个初始值，并通过迭代逼近方程的根。fzero函数的使用方法比较简单，但是对于复杂的方程可能需要多次尝试不同的初始值才能得到正确的结果。
f=fzero(@(h)(h+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./...(-10248.12+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./...(-10073.98+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/...(-10073.98+78.28+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/...(-10073.98+78.28*2+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/...(-10073.98+78.28*3+31649.66.*h).^2)+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).*...(sqrt(1+(-22143.12+11760-11760*r*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2)...-sqrt(1+(-23655.75+7.*9.8.*22.05+11760-11760*r*0.87-rou.*L.*9.8.*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+...2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2))./(rou.*9.8.*0.87)-H),0.9);
h=f;
phi0=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10248.12+31649.66.*h).^2));
phi1=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./(-10073.98+31649.66.*h).^2));
phi2=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28+31649.66.*h).^2));
phi3=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28*2+31649.66.*h).^2));
phi4=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.Z^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28*3+31649.66.*h).^2));
theta=asec(sqrt(1+(-23655.75+7.*9.8.*22.05-rou.*L.*9.8.*0.87+11760-11760*r*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2));
alpha=asec(sqrt(1+(-22143.12+11760-11760*r*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2));
x=sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta))./(rou.*9.8.*0.87).*log((sec(alpha)+tan(alpha))/(sec(theta)+tan(theta)))+sin(phi0)+sin(phi1)+sin(phi2)+sin(phi3)+sin(phi4);
%下面是判定锚链是否堆放在海床上，并进行此种情形的求解。
outFywrong=(-22143.12+11760-11760*r*0.87+31649.66*h)%输出错误情形下 Fy 的值
outGrou=rou.*9.8.*0.87.*L%输出错误情形下锚链的重力
if sign(rou.*9.8.*0.87.*L-(-22143.12+11760-11760*r*0.87+31649.66*h))==1
f=fzero(@(h)(h+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./...(-10248.12+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./...(-10073.98+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/...(-10073.98+78.28+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta))^2/...(-10073.98+78.28*2+31649.66.*h).^2)+1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/...(-10073.98+78.28*3+31649.66.*h).^2)+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).*...(sqrt(1+(-22143.12+11760-11760*r*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2)...-1)./(rou.*9.8.*0.87)-H),0.7);
h=f;
phi0=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10248.12+31649.66.*h).^2));
phi1=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2./(-10073.98+31649.66.*h).^2));
phi2=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28+31649.66.*h).^2));
phi3=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28*2+31649.66.*h).^2));
phi4=acos(1./sqrt(1+sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2/(-10073.98+78.28*3+31649.66.*h).^2));
theta=0;
alpha=asec(sqrt(1+(-22143.12+11760-11760*r*0.87+31649.66*h).^2/sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta)).^2));
x=sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta))./(rou.*9.8.*0.87).*log((sec(alpha)+tan(alpha)))+sin(phi0)+sin(phi1)+sin(phi2)+sin(phi3)+sin(phi4)+L-(31649.66*h-22143.12+11760-11760.*r.*0.87)./(rou.*9.8.*0.87);
end
outphi0=phi0.*180./pi%钢桶倾角
outphi1=phi1.*180./pi%钢管倾角
outphi2=phi2.*180./pi
outphi3=phi3.*180./pi
outphi4=phi4.*180./pi
h
outalpha=alpha.*180./pi
outtheta=theta.*180./pi
outx=x
outk=sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h).*(374.*2.*h.*V.^2).*cos(beta))./(rou.*9.8.*0.87)%锚链系数 k
outFy=(-22143.12+11760-11760*r*0.87+31649.66*h)
outFx=sqrt((2.5.*v.^2-1.25.*v.^2.*h).^2+(374.*2.*h.*V.^2).^2+2.*(2.5.*v.^2-1.25.*v.^2.*h))
放到MATLAB里运行，结果跟我说：
为啥呀？
