数学建模社区-数学中国

标题: 2016数学建模国赛A题程序（原创）作者cclplus [打印本页]
作者: 杨利霞    时间: 2019-4-10 10:54
标题: 2016数学建模国赛A题程序（原创）作者cclplus
2016数学建模国赛A题程序（原创）作者cclplus
3 f: a( h4 P6 L* m

2 L( n% t! m: v2 I
. z8 H2 b# O) j# e( [% E0 z& B! cclear all;) y0 B* L% Y" R. K; l
close all;
* n6 H# z7 z0 bclc2 b8 y2 d5 A0 H8 @( Y( o, I1 d
format long
) a, ?/ U7 q9 P4 Y6 h- }9 J' f- M; usyms h S Fw Ff Ff1 a b c d l L F depth n pl m x1 y1 y t distance n a1 b1;2 X5 g4 u$ ^' N0 |) Z$ Y* \5 m9 Q7 P8 K
F=[];
7 t1 w+ H7 f* x! C8 _theta=[]; - k3 [: j% g, i
v=24; %风速; C8 u; `1 L  f4 I+ o2 ?
l=105*10^(-3); %锚链每节链环的长度9 e" A* f$ \- |
L=22.05; %锚链的总长度
: G; Y" f- I5 ?* T, C) Nnum=0; %通过更改不在海床上的链节数得到一个最优解* E$ n. ?; T1 x8 O5 \; k0 o# r$ ]
num1=round(L/l);
# x" n/ z+ \- f& l% i* Q  mnum2=0;
7 h' C  v4 v3 l/ T3 m" `) Q) |' M6 Plin=0/180*pi; %第一个链节与水平方向的夹角
7 z5 e5 q3 H9 `  Y  k; nlin1=90/180*pi;: V. {  X/ `0 u1 z0 R
lin2=0;+ v, i! m) U2 ^, ^% A3 ^- F8 e% V
m2=1200; %重物球质量" k; R# C( B/ A+ [; n
pg=7.7*10^3; %重物球的密度（单位：kg/m^3)' C* z; Q6 T4 r- r! E
depth=20; %水深% e4 W0 E1 H$ R$ u0 k6 \8 K* C9 U
pl=7; %锚链单位长度的质量( o) c3 _  Y1 Q# p9 m" V! V
vh=0; %海水流速2 q! k2 a5 q9 {* C; e3 ~
g=9.8; %可通过改变此语句来修改重力加速度，单位为m/s^23 {% @( k( h$ o. s5 j6 N- ]" V
p=1.025*10^3; %海水密度
" T% U" M9 P7 W8 y" ?M=1000; %浮标质量
& L3 I. N, ~3 L  k( l- M' S6 L' d- E$ O8 ^m=10; %钢管质量8 |1 [/ N* r3 j! ^0 C2 \/ i
m1=100; %设备和钢桶总质量( a+ P9 J$ L4 f5 b* C# R" v, }
y=0;
: A# W% Y. ]$ B) F% O& X; R6 Yd=1;: F& H  s3 ]$ }/ @% {8 H2 W) n
j1=0;5 K- y% _% Z+ |7 f' `, K" Y
j2=0;
8 ~* g. G3 U3 t- o: @( qwhile(abs(y-d)>0.005)%在这里选择所需要的精度，
  ~& z' `4 w" V" D9 _if (y>d)&&(num<round(L/l))
  ]! D; [" w. knum1=num;. {# F' K; q/ X! x, ^3 z$ O/ A
num=round((num1+num2)/2);
5 k+ x9 m; \4 J5 D! S% ?; Melseif (y<d)&&(num<round(L/l));
/ M+ H5 Q& s/ K; B& N% ^num2=num;+ F( B3 w( L( T$ X( i
num=round((num1+num2)/2);
$ _/ F. H$ p2 b1 [# Melseif (y<d)&&(num==round(L/l)); R2 `) N" v7 G
lin2=lin;! q' g+ i& v: F4 g* c
lin=(lin1+lin2)/2;+ m- `9 h% F' t8 ^" j. a( j
elseif(y>d)&&(num==round(L/l)); F1 d3 U" |. D- ^9 k: @7 M
lin1=lin;
" G7 Z2 C; x- i: Y2 ?+ llin=(lin1+lin2)/2;
" T2 y0 E/ ~3 T+ G0 J% R" B! v2 j( n; jend
  u6 y+ z: l4 O" x0 B. z%钢桶受到的浮力! ^2 ~, K: k: j4 `' p
Ff1=p*g*pi*(0.3/2)^2;! _7 m# v' o6 ]& z4 ~
%钢管收到的浮力
* G1 G. q& |& @6 A* MFf2=p*g*pi*(0.05/2)^2;
5 P6 A" r5 e/ ~: D5 R%重物球所受浮力
7 H  v# y' M9 N5 ^# u3 ^Ffg=p*g*m2/pg;  R% M$ U5 O1 E3 Q& Y3 x2 J
%重物球所受海水水流力0 v1 t* q. f/ i. S& S
Fhg=374*pi*((m2/pg/3/4)^(1/3))^2*vh^2;
" A0 _5 {/ K: J; s! z9 r; J%风对浮标受力面的投影面积/ p( [9 u  N' G+ p
S=2*(2-h);' w- b+ ~( X* X& v: x/ J7 A% d
%风对浮标产生的力+ U5 m! U; j  r( h8 W5 t
Fw=0.625*S*v^2;* ]# J9 x2 h- [
%浮标在水中的体积2 ?/ Z" C" z6 I6 z, F* t; L: v% q
V=pi*(2/2)^2*h;; X$ z: f, ]3 T2 g8 p
%浮标所受到的浮力) J& q3 ~2 x, F( k  i3 x: v, C
Ff=p*g*V;( w3 i, s, @3 r' j, C
%浮标受到海水的近似水流力' @5 E* C  D5 m# k5 f
Fb=374*2*h*vh^2;
- c8 o+ P; q1 k2 w* r! M- Q%钢桶受到海水的近似水流力
0 G5 ~0 d# \! p* u) g7 M' u! TFs1=374*0.3*vh^2;
, i. X( ?  q- R8 j! G%钢管受到海水的水流力的近似值
. i2 `, }; t0 ?: Q7 ]1 zFs=374*0.05*vh^2;+ P: X4 u& M% r4 y3 c, n" M
%浮标浸没水中的高度  ~4 O+ y; u0 O4 Z7 u
if num==round(L/l)( C. v* U& d( i" e1 A
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));( M' ^: [( g, y2 q3 W3 [
else % w9 L/ @+ T8 a; n# l/ o7 B
h=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+num*pl*l*g)/(p*g*pi);" u  X" x) ]# G5 P- R+ [+ [
end8 U/ w: l0 n+ d  o
a=Fw+Fb;
8 r$ Z7 I) y7 W' Wb=-M*g+Ff+(Fw+Fb)*tan(lin);
8 \" Q; y( G1 S+ Z0 v% r/ Z% t! |if j1==05 W6 y% ?/ z. W' m" \# ?
a=eval(a);# f5 f$ J* }3 R4 V" q0 C
b=eval(b);8 c* `* y/ n$ E& `
else$ ?! l9 A% @0 Z7 Y
end- g) h7 i0 z% ^; L6 h; P
F(1)=sqrt(a^2+b^2);
3 _4 ]1 w( I: p8 C1 otheta(1)=atan(b/a);) O/ W3 j. `5 m( ?
n=0;
& X1 ^( D  n$ C  f% G3 h. k. vfor i=1:4; i; ~3 S* }8 P1 ]1 K0 s5 q
%钢管受到海水的水流力( N8 T9 ^" U% b
Fh(i)=374*0.05*sin(theta(i));
! |: v: X( A3 U6 O0 A8 Ln=n+Fh(i);
% ?( B& ]' t! Q7 ?a=Fw+Fb+n;
9 h8 \8 M7 A4 i; uif j1==0
9 f# e0 d& g: k0 Q9 ta=eval(a);: E; w3 v- Y% s6 |" [6 T
else- _  ?3 I2 Z3 P2 I( w5 F3 r
end! v+ Y/ C7 G: I# a4 n
b=F(i)*sin(theta(i))+p*g*pi*(50*10^(-3)/2)^2-m*g;% T1 F5 ~7 c* k) M6 ~
F(i+1)=sqrt(a^2+b^2);2 S% F0 q8 h/ N6 Q( U/ m
theta(i+1)=atan(b/a);  z$ V- @. y, D; ?
end
0 i- c1 i* u. d- I! W. ?c=0;
- I$ h8 @6 I7 ~$ e" Afor i=1:5
8 Z6 a) C; J- _) d( v4 R7 f7 Sc=c+sin(theta(i));
9 g9 r. b0 W: x5 M; Z: ~8 U& Z( Hend
$ v. d) f$ c3 m' c- f  y& S; ~! |" md=depth-c-h;7 b! I$ L/ ~" _* X3 ^9 B" \( P
y1=lin;
8 Z2 q0 \8 z0 ldistance=0;
% x) X0 G! b+ c2 Q0 @if num==round(L/l)
. s! [( x3 v5 k& M! j$ ]y=l*sin(y1);, ~& S! G4 Q+ M# r7 V: U
x1=Fw/sqrt(1-(sin(y1))^2);
! X4 u: s7 `- \" Q) X/ afor i=1:num-13 g  x0 T- m+ W5 Y! T
m=(x1*sin(y1)+i*pl*l*g)/sqrt((x1*sin(y1)+i*pl*l*g)^2+Fw^2);- v1 [- g" u; b: U6 H3 x
m=m*l;
6 G6 z" A* F; u" }& \. E% Vy=y+m;* z9 }. [1 q/ f. v, l- o- F# ]
n=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;$ I/ T# R+ F+ r$ r
if j1==0
: h; {: R' C2 l5 ~n=eval(n);0 k# _- k/ e: l) _4 w# T( {
else6 f  b$ Z1 ~3 a: y  g! n
end" c' j4 i; r% z7 c+ U, H
distance=distance+n;5 j- ^& q" o3 l% J! G
if j1==01 Y( V6 W3 C6 O3 J* r; C. U/ d
y=eval(y);3 K+ c( \% ~6 p+ q4 p
else
  l. I3 c' I6 }/ H6 aend
7 s/ V1 N/ _+ ~$ L) Rend' C" i* ~4 ?! \1 S) i; t; L
else
. G3 A' j+ I. |$ e3 E0 e. |y=y1*l;
( z! M7 J% q3 I+ v* udistance=(round(L/l)-num)*l;
, G, N5 F2 w; V( i1 a& qfor i=1:num7 x/ o( }0 d: s1 c3 t: E2 z
x1=Fw/sqrt(1-(sin(y1))^2);0 |9 {- r, V  n
m=(x1*y1+i*pl*l*g)/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
& [) b" n7 l! Y% t" ey=y+m;
9 a& @1 a6 B4 z0 yn=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
) }7 q# s$ B. |8 xif j1==0: ~% i1 U; o% m9 M! F  g
y=eval(y);8 a9 j' l- |8 a) t( p
n=eval(n);1 }! D& h- C8 z/ j6 Y' o3 o
else' K0 t) b4 [/ g* _6 v" e
end
( z/ b, W' G7 h) _$ Odistance=distance+n;+ g5 A' U( t. f+ `- _) ?
end
0 ^3 Y  @# `0 M) aend4 L# O, l6 N4 |$ w
m=0;8 b# T& K2 m  U2 c5 k1 u2 T& E
j1=1;
& t- o+ N) s  G& _; a- }8 bj2=j2+1;
7 u0 r9 m) w1 M7 R! v8 j$ I1 Aend+ P/ I: |& M& F! F
%钢桶受到的浮力/ V# W3 J; N5 v8 s/ ]1 U
Ff1=p*g*pi*(0.3/2)^2;
6 Y  V" ^2 {# ?% c8 T$ x- C. Y% ^%钢管收到的浮力5 a* `- T7 I3 ^: h7 {' q. O1 P
Ff2=p*g*pi*(0.05/2)^2;7 T# A. `4 i+ {4 E' t7 N
%重物球所受浮力
3 z! F! m4 b" C3 j$ g8 h: yFfg=p*g*m2/pg;) {! ]0 O1 x4 _+ [5 F
%重物球所受海水水流力
9 R5 B* i! b( n! l' QFhg=374*pi*((m2/pg/3/4)^(1/3))^2*vh^2;: ?9 W8 l$ ?. I
%风对浮标受力面的投影面积
! J, F/ s4 y* V, r1 h/ xS=2*(2-h);
$ s' C9 f- k; I+ w7 U7 T  n%风对浮标产生的力
( d* g' L+ c8 k  D/ V/ LFw=0.625*S*v^2;: t3 d2 j/ C: @' \+ ]( x! M
%浮标在水中的体积
5 Q2 i/ L) ~. [4 x. e( y+ wV=pi*(2/2)^2*h;9 F* `% Y4 h# A3 |# M; ^
%浮标所受到的浮力
+ F  r* p: q! |: E. RFf=p*g*V;/ z$ r# R8 k# s0 z0 v; w8 M* I
%浮标受到海水的近似水流力
9 {0 F  ^& ~% T) j6 M' `( j/ ]& qFb=374*2*h*vh^2;
9 A+ r; C* R. G: N%钢桶受到海水的近似水流力
# m; K- D$ ?# G" q7 OFs1=374*0.3*vh^2;
$ {) x" t3 A8 ?$ w" T  y! M4 S  x%钢管受到海水的水流力的近似值/ G7 W& T% `1 O1 g% r* Q+ F( I
Fs=374*0.05*vh^2;
' j: [. ~. w. Z0 J/ K9 p%浮标浸没水中的高度
2 o" H) d% ]3 ?) X( xif num==round(L/l)
# v; u$ Q: q* x5 ]3 `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));
8 U7 h: s7 H, [0 D# P9 b& ?# selse
3 P: L3 D7 v, Y$ }# Rh=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+num*pl*l*g)/(p*g*pi);0 [* P+ |% u% @- D4 Q
end* M& e' Y7 f4 M8 U* g. ?) ]  }& P% E
a=Fw+Fb;6 M  s' p6 a( R1 b. I6 G! C
b=-M*g+Ff+(Fw+Fb)*tan(lin);7 e- Z6 r6 ]9 p0 C% A
F(1)=sqrt(a^2+b^2);
* X  _1 b) s! T+ E6 z$ o4 t& c8 H/ h5 Xtheta(1)=atan(b/a);
0 z7 O; Y3 x; m* @8 }) |1 ~1 q5 S/ ln=0;
) D* E5 T" Z8 ]3 |  Zfor i=1:4+ e9 m2 ~! M7 k3 y* ]5 A* p
%钢管受到海水的水流力& q! P, C" K! A) q/ b% s
Fh(i)=374*0.05*sin(theta(i));
$ D( U% H% i' I: R; ]# kn=n+Fh(i);
  G# ~- o+ h1 k. ^* S1 o  ia=Fw+Fb+n;
& X: Q- O$ V! |. z, gb=F(i)*sin(theta(i))+p*g*pi*(50*10^(-3)/2)^2-m*g;- b& [$ z! z" w! h3 g( i/ H
F(i+1)=sqrt(a^2+b^2);9 t: S: f* P) r" C/ X
theta(i+1)=atan(b/a);
' k% F, [' o1 N) G; k/ K1 Wend: ?; e# O, ^: l
disp('输出钢管和钢桶的倾斜角度（角度制）'), @& V; D4 V' j! {
th=90-theta*180/pi& @$ b- F9 L; ^% B8 ?
m=85*pi/180;
, U! q- h$ I5 Q1 O, H4 O; t8 Tif theta(5)>m
+ M5 k' Z+ B. k# i: }/ A4 mdisp('钢桶的倾斜角足够小，测量准确')
& B$ U* b+ |5 C0 }else 8 Y( z8 |) h) S) s  s7 M
disp('钢桶的倾斜角过大')
/ W. e; q1 _! c% r; e' aend+ M2 U1 r+ F# h; X' i
c=0;" i  K3 D8 s$ o; P
for i=1:5
+ L, `3 ]" n2 ]& o6 vc=c+sin(theta(i));- H, f& n3 _' P$ o9 u
end
. W6 f/ A2 H0 q/ D; W/ fd=depth-c-h;' W, ~* g" z' ^& M0 D
y1=lin;$ c, V6 v6 J* L) P4 x
distance=0;6 p# b9 |. R5 c! C! S
if num==round(L/l)
! c. H7 u. b+ B. {1 p- j$ ?y=l*sin(y1);3 t( g) ^  s% |8 B' E4 i* _. U3 F
x1=Fw/sqrt(1-(sin(y1))^2);
' w$ |* l' ~& |3 k1 v" Wfor i=1:num-1
$ V$ B/ G0 g% f+ ~8 um=(x1*sin(y1)+i*pl*l*g)/sqrt((x1*sin(y1)+i*pl*l*g)^2+Fw^2);
6 A. V6 {4 H) Xm=m*l;5 _% d/ d! ?' I# |+ O
y=y+m;. u+ b$ L3 B6 A. U
n=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
( p. \# X  P0 Zdistance=distance+n;7 p9 M* O4 {, M2 ~3 z
plot(distance,y,'o')
0 ]0 P8 Q7 e6 j3 A) U; f$ mhold on
3 o1 t4 U8 S" [end
& K! {8 p! d4 d+ Lelse3 u7 y* }& h6 p3 U+ n6 g
y=y1*l;+ R, v. q% h. i  |
for i=1:round(L/l)-num5 L4 t( j9 {. e2 c5 O7 O
distance=i*l;1 E6 \" U- z0 z2 o! Q$ |4 O  A
y=0;* J* V# Q3 |2 F
plot(distance,y,'o')' I) ~! m! z- Z
hold on
4 i# U; g! N  H+ Pgrid on8 |8 p6 _6 z4 i3 V0 q2 K  i
end+ A! l6 V( z; ?1 v" Y! s* @! R
for i=1:num
& b, n8 t1 _$ Q% \0 Mx1=Fw/sqrt(1-(sin(y1))^2);/ i" w0 o) K% W7 s  _6 z
m=(x1*y1+i*pl*l*g)/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;6 B( P: Z, V+ e% W8 K! J3 m- J
y=y+m;
, D. B$ U/ F7 ?' P5 Gn=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
! L- E7 I0 u  C+ {if j1==0# p% U' O" T6 a4 c) S! \2 a$ d
y=eval(y);
2 s) R% k. `& m7 Pn=eval(n);) i, x" z% G& t; |
else
7 U% R9 M' g! h! Q! ~0 ?end. F( K% g, p. d  a7 v8 P
distance=distance+n;. Y" D4 e0 U5 y* _* W5 }
plot(distance,y,'o')
9 P/ b" z: P5 Q' b- V# \5 Ehold on
/ |4 H! G) F# H, wend% Q+ d; ~( K: g9 d
end
! Z5 Y$ a1 z2 |  C& cm=0;; r2 Y( r( b. s8 q' @8 [
for i=1:5; y4 \6 q/ R5 J; X+ f; H
m=m+cos(theta(i));
# l/ m$ x$ d% `/ q7 Iend6 `; a/ E& g1 c* |; J
%浮标的运动半径
: g# z8 q- n4 O" ]5 _disp('输出浮标的运动半径')
5 B: k' K- z3 O- A. @7 }ans=distance+m9 t; c) W! ~. B' F4 m9 Z
( D0 Q! [8 q$ W/ s2 L1 x

3 [/ b6 X: t) _# r/ [! b, z
: R0 M$ O6 c5 R% r3 U2 V2 j0 X4 L3 q7 Z0 \* h' v

$ A$ s) o' `. D# m' z2 {
6 [1 x( o  ]# s& S
: ?7 u, y' f8 o* x1 d9 {

2018全国数学建模总结.docx

17.26 KB, 下载次数: 0, 下载积分: 体力 -2 点
作者: 571334077    时间: 2019-4-10 19:25
2333333333333333333333# _  y. \! u. @% {* r




欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) Powered by Discuz! X2.5