2016数学建模国赛A题程序（原创）作者cclplus

杨利霞

2016数学建模国赛A题程序（原创）作者cclplus

! W. Y& X  u& s9 K, c
* r  M0 ]0 \: T2 u
0 @4 x6 f5 F. k9 B; p0 i7 P% w% Kclear all;
! C3 F1 }, _" i3 _. z: d: ]close all;
clc
7 f* {+ h  u; ?( Z, b2 Sformat long
syms h S Fw Ff Ff1 a b c d l L F depth n pl m x1 y1 y t distance n a1 b1;
F=[];
) @7 x4 T* P; k( a( Q* ?: ^3 `  @! Dtheta=[];
$ |* k6 h2 M) Y9 \& Cv=24; %风速
l=105*10^(-3); %锚链每节链环的长度
+ X7 G+ L5 O; N+ q% \, O7 G* ]L=22.05; %锚链的总长度
num=0; %通过更改不在海床上的链节数得到一个最优解
! D3 u" W6 W+ \$ }" d) s" lnum1=round(L/l);
; E# y% D5 [7 I  N  ~1 e3 ^/ lnum2=0;
lin=0/180*pi; %第一个链节与水平方向的夹角
9 R+ u" d8 [" i* `lin1=90/180*pi;
  r$ V  ?; s2 i* t: `lin2=0;
( H1 f/ ?5 N9 Y# P& a; sm2=1200; %重物球质量
* _) W( A, z1 e* e4 @% Lpg=7.7*10^3; %重物球的密度（单位：kg/m^3)
( [' B, f2 k% s  @0 R/ e1 adepth=20; %水深
% d- l) U/ }8 ]2 V! M1 Apl=7; %锚链单位长度的质量
vh=0; %海水流速
g=9.8; %可通过改变此语句来修改重力加速度，单位为m/s^2
p=1.025*10^3; %海水密度
6 ]- o) V- ?1 B  `. Z: }+ `M=1000; %浮标质量
: i$ t4 J; _' j, l" G7 |' Tm=10; %钢管质量
* `5 `" y& q: y  e, ^m1=100; %设备和钢桶总质量
y=0;
d=1;
. `  x% x7 Q+ r" V. W( d2 t+ Ej1=0;
j2=0;
while(abs(y-d)>0.005)%在这里选择所需要的精度，
if (y>d)&&(num<round(L/l))
6 Y: p" G. V/ }- Hnum1=num;
& V- ?8 ]6 l' p; ]2 y; l# unum=round((num1+num2)/2);
6 d- R  |& `7 E* U" C; {elseif (y<d)&&(num<round(L/l));
& {8 O! z4 d8 ~5 |5 q* n3 \0 ^num2=num;
# U9 N$ Q. P2 a( r  hnum=round((num1+num2)/2);
elseif (y<d)&&(num==round(L/l))
  X0 c  G3 I0 y1 r- l% m9 vlin2=lin;
lin=(lin1+lin2)/2;
elseif(y>d)&&(num==round(L/l))
lin1=lin;
lin=(lin1+lin2)/2;
$ y/ d5 Z6 ~/ Uend
  X% I. R! v! b* X%钢桶受到的浮力
Ff1=p*g*pi*(0.3/2)^2;
- e  W# z5 a0 |/ f% Z2 p%钢管收到的浮力
' }0 w! q9 F- @: ?7 i+ w. l- |% Y( cFf2=p*g*pi*(0.05/2)^2;
4 q9 z# x$ O0 a( C%重物球所受浮力
- c% i8 d$ v* ?/ P' O) mFfg=p*g*m2/pg;
%重物球所受海水水流力
' [- U% Y* f$ o; F0 vFhg=374*pi*((m2/pg/3/4)^(1/3))^2*vh^2;
%风对浮标受力面的投影面积
S=2*(2-h);
%风对浮标产生的力
Fw=0.625*S*v^2;
%浮标在水中的体积
: d3 T2 ]) ], S0 z, }V=pi*(2/2)^2*h;
%浮标所受到的浮力
Ff=p*g*V;
%浮标受到海水的近似水流力
4 S5 R5 D* s; Y/ ]+ e4 CFb=374*2*h*vh^2;
%钢桶受到海水的近似水流力
Fs1=374*0.3*vh^2;
%钢管受到海水的水流力的近似值
Fs=374*0.05*vh^2;
6 f6 }7 j2 d% C7 L0 e& P%浮标浸没水中的高度
if num==round(L/l)
h=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+pl*L*g+(Fhg+4*Fs+Fs1)*tan(lin))/(p*g*pi-(1.25*v^2+374*vh^2)*tan(lin));
else
h=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+num*pl*l*g)/(p*g*pi);
end
: I) I: `- ~9 |7 N& ^a=Fw+Fb;
b=-M*g+Ff+(Fw+Fb)*tan(lin);
2 {& Q9 J/ W$ m8 bif j1==0
* W4 ]* @' Y5 \' m0 Ba=eval(a);
b=eval(b);
/ W; r1 q- {2 ^* d* J# V/ x+ relse
) w/ w6 Z2 {/ g$ a/ ]end
F(1)=sqrt(a^2+b^2);
theta(1)=atan(b/a);
1 a$ [7 Y) a. o# C4 V$ An=0;
+ j! l+ o9 t9 s! T6 b  ofor i=1:4
%钢管受到海水的水流力
7 n" v0 y( J( u+ O# YFh(i)=374*0.05*sin(theta(i));
+ O9 U3 ]# V0 c8 Y) T7 g2 _n=n+Fh(i);
. K5 y+ s4 q( o, ]  E. v! k) Ca=Fw+Fb+n;
7 K- g* I/ ]! H$ wif j1==0
a=eval(a);
# u2 H; J. }% A0 b. s  N# K/ jelse
end
8 V  N5 ?! Z; y6 U6 I0 E4 _- db=F(i)*sin(theta(i))+p*g*pi*(50*10^(-3)/2)^2-m*g;
F(i+1)=sqrt(a^2+b^2);
theta(i+1)=atan(b/a);
! c9 a; U& x2 hend
6 W$ u6 T4 P6 [- i8 ]/ bc=0;
for i=1:5
: W" c3 S$ ?) b, I/ T$ q8 ^c=c+sin(theta(i));
end
% ?. q9 c) ?0 O9 O8 Dd=depth-c-h;
" l$ `, `2 D/ \- O9 Ey1=lin;
distance=0;
if num==round(L/l)
y=l*sin(y1);
x1=Fw/sqrt(1-(sin(y1))^2);
for i=1:num-1
m=(x1*sin(y1)+i*pl*l*g)/sqrt((x1*sin(y1)+i*pl*l*g)^2+Fw^2);
m=m*l;
$ C, V2 o* X4 Y' Y9 Ny=y+m;
7 \( a5 Y& j% F; x% Sn=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
3 I/ p; U) ?9 f- Z6 z+ Qif j1==0
) A. S2 r* V8 r6 Sn=eval(n);
0 H. A( r1 k6 U) x1 y2 d( Zelse
end
distance=distance+n;
if j1==0
  {" C9 s$ B/ r$ \1 X& T$ Ay=eval(y);
% d7 c6 ?  O/ r* |; A( y- l0 m% Uelse
end
end
, L. Q8 _6 s- _* o1 b8 L3 helse
1 f. T0 u$ S1 d' j7 my=y1*l;
distance=(round(L/l)-num)*l;
' l) v8 n4 K5 O, O3 V/ F/ `for i=1:num
4 a; x2 D1 b6 G* }/ D# q4 xx1=Fw/sqrt(1-(sin(y1))^2);
% w* U9 C  _8 \7 J; y  |2 um=(x1*y1+i*pl*l*g)/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
0 e% a  m/ `4 h# f, t# W9 Y' s) ~y=y+m;
n=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
( _; n4 f5 q3 A3 O4 e/ ]if j1==0
: b9 R* ^2 \( C: L7 L+ fy=eval(y);
8 Z% J- R$ J3 T. {" Sn=eval(n);
else
! `: T- ~5 H+ dend
- M6 t4 m" C6 L0 J% Q" F2 m) r+ edistance=distance+n;
/ l3 T% ^) J2 `! ]end
* z3 _0 m, U1 g+ E' P9 r% l" B7 Rend
! r/ @7 j2 g+ f, Wm=0;
9 y! X4 q* y' S( [j1=1;
j2=j2+1;
end
7 g# |; \' W7 x; a%钢桶受到的浮力
$ Q( T, e3 M! c7 @8 T! L% sFf1=p*g*pi*(0.3/2)^2;
' b( E9 z) F2 r3 c7 n) U- R%钢管收到的浮力
Ff2=p*g*pi*(0.05/2)^2;
( L8 [$ w0 [# \2 x8 A& S%重物球所受浮力
Ffg=p*g*m2/pg;
%重物球所受海水水流力
, w& G7 R- ]* J1 ^  ~Fhg=374*pi*((m2/pg/3/4)^(1/3))^2*vh^2;
%风对浮标受力面的投影面积
S=2*(2-h);
/ x" A7 u8 j5 c9 ]2 m1 k- C4 h) i+ c  M& i%风对浮标产生的力
Fw=0.625*S*v^2;
3 r* V$ o0 s' M7 j# M+ K; k( `%浮标在水中的体积
( E4 m- U3 @# A4 oV=pi*(2/2)^2*h;
%浮标所受到的浮力
/ X$ z3 R& X0 a, q: d5 ?8 Z+ [4 @! sFf=p*g*V;
%浮标受到海水的近似水流力
" f1 ^% `& Y4 J. E( ^6 ?Fb=374*2*h*vh^2;
%钢桶受到海水的近似水流力
; z; W7 Z) s/ y; K; b, eFs1=374*0.3*vh^2;
: N% Z3 ?3 E6 E! M: h%钢管受到海水的水流力的近似值
Fs=374*0.05*vh^2;
/ G& S& \' V' m" _/ h2 j! M%浮标浸没水中的高度
6 U" e% o2 E- w" bif num==round(L/l)
h=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+pl*L*g+(Fhg+4*Fs+Fs1)*tan(lin))/(p*g*pi-(1.25*v^2+374*vh^2)*tan(lin));
# u- |+ [. r! k( k1 gelse
h=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+num*pl*l*g)/(p*g*pi);
end
1 P9 x; I" L/ v# ia=Fw+Fb;
) Y2 z2 w+ V- A, Y% Xb=-M*g+Ff+(Fw+Fb)*tan(lin);
8 M3 \1 e9 g! O# z& tF(1)=sqrt(a^2+b^2);
theta(1)=atan(b/a);
n=0;
, n7 i: @+ m( N, S" Cfor i=1:4
%钢管受到海水的水流力
; J9 l7 W4 L  [3 X1 EFh(i)=374*0.05*sin(theta(i));
n=n+Fh(i);
) ~4 T* k$ W) ~% p. a$ @7 d' ~a=Fw+Fb+n;
b=F(i)*sin(theta(i))+p*g*pi*(50*10^(-3)/2)^2-m*g;
- c4 B6 u. K" Z) q: T  _F(i+1)=sqrt(a^2+b^2);
' F$ L% R# P9 ~# ttheta(i+1)=atan(b/a);
end
disp('输出钢管和钢桶的倾斜角度（角度制）')
4 @& t9 C+ h% q) A7 l  Yth=90-theta*180/pi
2 D* ], |1 B4 ^, jm=85*pi/180;
& u; ^. w& g! X" Zif theta(5)>m
disp('钢桶的倾斜角足够小，测量准确')
else
. w! f! `. S) Z) Vdisp('钢桶的倾斜角过大')
end
% f' g( V1 T- m& R% r) H: Bc=0;
% u! v0 g* J3 I* v! v' V" T# yfor i=1:5
c=c+sin(theta(i));
4 m" ]- t( x5 k; x2 rend
& ^! w2 S& P$ T  d9 f/ [d=depth-c-h;
y1=lin;
1 x4 z# c9 D( F% E5 o- d+ Ldistance=0;
# r, V3 U* \. Hif num==round(L/l)
4 p% W2 |; _# d( V- Q9 {y=l*sin(y1);
  u7 M3 g3 P% Q  k( S' }x1=Fw/sqrt(1-(sin(y1))^2);
/ |: ?  J0 S$ B; c1 Gfor i=1:num-1
8 @' z* R3 `1 q0 i- o. @, g7 |6 tm=(x1*sin(y1)+i*pl*l*g)/sqrt((x1*sin(y1)+i*pl*l*g)^2+Fw^2);
- O) d: K4 ]6 ]- T$ ~3 }& G! vm=m*l;
( T9 H2 |3 r: s- q' ?  q5 Qy=y+m;
5 v+ y9 D0 o8 c  x: O, Vn=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
- x2 j$ J$ w/ ?: V" K$ g2 e. adistance=distance+n;
% }$ v/ i* z. {! o- O9 kplot(distance,y,'o')
8 n3 @! Z: k0 X/ _. ]5 jhold on
7 L0 y- T( {: G. I$ M2 J2 rend
6 M3 X8 r: {" N5 O; s* G1 ^else
$ M- [0 G2 {; Y3 O/ u# hy=y1*l;
4 }7 V( U: J1 J( C7 ofor i=1:round(L/l)-num
  I0 R  l3 a/ Cdistance=i*l;
; L$ ]2 j+ A$ C) O: Jy=0;
+ Z* x4 M# D1 L/ i/ vplot(distance,y,'o')
$ O& {" S9 z9 J, ^. S. m/ bhold on
8 X, }3 u, n* [6 zgrid on
- \- d5 L$ w& M; Xend
( M- i- b; ]7 _2 \2 |) Z0 f. ^0 Rfor i=1:num
" z' t: e8 e" W' Xx1=Fw/sqrt(1-(sin(y1))^2);
$ c% g9 ^4 S' S$ T; om=(x1*y1+i*pl*l*g)/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
/ a+ E0 [4 e4 |9 M6 O; n+ g% by=y+m;
7 s- W0 r1 c7 c% `! s& N! A. f# ln=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
. v7 H- t. @6 I6 u4 ?. [if j1==0
y=eval(y);
/ ]8 J0 s! M' {4 c8 Hn=eval(n);
5 W/ ?0 i3 a, k' s: r5 relse
7 k% n$ l+ C! g5 |end
distance=distance+n;
plot(distance,y,'o')
2 h$ m& J/ G' z! T2 |& N' Rhold on
" ]2 H7 t- h% K8 Tend
6 e. C  q% h4 wend
& |# |$ ^% a6 m, s- S, Lm=0;
% ]# v, ]6 u9 R  V8 Q  r' ofor i=1:5
m=m+cos(theta(i));
end
9 X4 @' }5 }- o" `%浮标的运动半径
disp('输出浮标的运动半径')
ans=distance+m
. I) l: o  S) L! N, v# E


: Y3 C! `8 r& J( n4 f, V# ]% ~1 {, T
! s8 f0 S: s& [# ]' b6 a4 z# J


! O$ e$ }8 G- W6 f( d5 S' d7 e# S/ M

571334077

2333333333333333333333