Matlab数学建模：连续模型

杨利霞

9 ?1 A2 F/ a5 a, V5 WMatlab数学建模：连续模型一、学习目标

（1）了解连续模型及解题步骤。

（2）掌握用 dsolve 求解常见的微分方程解析解。

（3）掌握用 ODE 家族的求解器求解数值解。

（4）掌握使用专用的求解器求解。

二、实例演练（1）谈谈连续模型在数学建模中的应用。

       连续模型是指模型是连续函数的一类模型总称，具体建模方法主要是微分方程建模。微分方程建模是数学建模的重要方法，因为许多实际问题的数学描述将导致求解微分方程的定解问题。把形形色色的实际问题化成微分方程的定解问题，大体上可以按以下几步：

    1. 根据实际要求确定要研究的量（自变量、未知函数、必要的参数等）并确定坐标系。

    2. 找出这些量所满足的基本规律（物理的、几何的、化学的或生物学的等等）。

    3. 运用这些规律列出方程和定解条件。

        MATLAB 在微分模型建模过程中的主要作用是求解微分方程的解析解， 将微分方程转化为一般的函数形式。 另外， 微分方程建模， 一定要做数值模拟， 即根据方程的表达形式， 给出变量间关系的图形， 做数值模拟也需要用 MATLAB 来实现。

微分方程的形式多样，微分方程的求解也是根据不同的形式采用不同的方法， 在建模比赛中， 常用的方法有三种：

    1. 用 dsolve 求解常见的微分方程解析解

    2. 用 ODE 家族的求解器求解数值解

    3. 使用专用的求解器求解

（2）常规微分方程的求解

微分方程在 MATLAB 中固定的表达方式， 这些基本的表达方式如下表所示：

对于通常的微分方程， 一般需要先求解析解， 那么 dsolve  是首先考虑的求解器，因为dsolve 能够求解解析解，其具体的用法如下：

[实例]求微分方程 xy'+y-e^x=0 在初始条件 y(1)=2e 下的特解，并画出解函数的图形．

求解本问题的 Matlab 程序为：

• 
  2 e; m. i! u' W7 J: o. Z
  syms x y
  
  
• 
  
  - b$ Z4 Q) K* P, c! l$ Uy=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x')
  
  
• 
  
  ezplot(y)
  + z. \/ L6 H1 c, N2 l4 [/ {
  2 C* p4 C. e( X2 y2 w( N
  

微分方程的特解为：y=1/x*exp(x)+1/x*exp (1) (Matlab格式)，即 y=(e+e^x)/x ，此函数的图形如下图：

（2）ODE 家族求解器

如果微分方程的解析形式求解不出来， 那么退而求其次的办法是求解数值解， 那么这个时候就需要用 ODE 家族的求解器求解微分方程的数值解啦。

因为没有一种算法可以有效地解决所有的 ODE 问题，为此，MATLAB 提供了多种求解器，对于不同的 ODE 问题，采用不同的 Solver。MATLAB 中常用的微分方程数值解的求解器及特点如下表所示。

要特别提醒的是：ode23、ode45 是极其常用的用来求解非刚性标准形式一阶常微分方程(组)初值问题解的 Matlab 的常用程序，其中：

ode23 采用龙格-库塔2 阶算法，用3 阶公式作误差估计来调节步长，具有低等的精度．

ode45 则采用龙格-库塔4 阶算法，用5 阶公式作误差估计来调节步长，具有中等的精度．

[实例]导弹追踪问题

设位于坐标原点的甲舰向位于 x 轴上点 A(1, 0) 处的乙舰发射导弹，导弹头始终对准乙舰。如果乙舰以最大的速度 v0 (是常数)沿平行于 y 轴的直线行驶，导弹的速度是 5*v0，求导弹运行的曲线方程，以及乙舰行驶多远时，导弹将它击中？

记导弹的速度为 w，乙舰的速率恒为 v0。设时刻 t 乙舰的坐标为 (X(t),Y(t))，导弹的坐标为 (x(t),y(t)) 。当零时刻，(X(0),Y(0))=(1,0)，(x(0),y(0))=(0,0)，建立微分方程模型：

因乙舰以速度 v0 沿直线 x=1 运动，设 v0=1，w=5，X=1，Y=t，因此导弹运动轨迹的参数方程为：

MATLAB 求解数值解程序如下：

(1) 定义方程的函数形式：

• 
  7 r9 `8 z# R; \( k& U/ t
  function dy=eq2(t,y)
  7 O2 U  n  h7 C: S) ^8 j
  
• 
  
  dy=zeros(2,1);
  
  ( F7 L0 ^1 g$ o
• 
  % h. q9 {2 m8 i/ M$ s2 s8 g
  dy(1)=5*(1-y(1))/sqrt((1-y(1))^2+(t-y(2))^2);
  
  
• 
  
  dy(2)=5*(t-y(2))/sqrt((1-y(1))^2+(t-y(2))^2);
  
  
  , U$ q. {1 l9 Q$ A1 [

2) 求解微分方程的数值解

• 
  2 p) E3 p0 O7 |  e; N% G% C
  t0=0,tf=0.21;
  
  
• 
  
  [t,y]=ode45('eq2',[t0 tf],[0 0]);
  6 V7 D4 f$ p" f2 W+ ~0 y
  ) l9 w) C0 |2 ]8 w3 U
• 
  
  / r  @1 c  Y, |& P, S6 R& LX=1;Y=0:0.001:0.21;plot(X,Y,'-')
  ! O4 \. `# Y& z: c; r& G, U
  
• 
  # \) J. K/ b  h4 e& J2 W' l, C
  ! f+ |# \' j! Y( I2 Xplot(y(:,1),y(:,2),'*'),hold on
  
  
• 
  
  x=0:0.01:1; y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24;
  
  
• 
  
  . b7 }$ f3 F! n' rplot(x,y,'r')
  
  
• 
  : W1 A1 v& R' v2 ?1 T8 Q
  - s. v0 s0 w. H- F4 |
  
  . r" |/ f$ j5 _3 j2 h
  

（3）专用求解器

对于复杂的微分方程模型的求解， 可以借助 MATLAB 偏微分方程工具箱中的专用求解器。以下将以一个实例来看看如何借助偏微分方程工具箱来实现一个微分方程的求解与数值仿真。

所研究的对象是一个二阶波的方程：

这个时候要查看一下 MALTAB 中哪个函数能求解相类似的方程， solvepde 可以求解的方程形式为：

' G( c/ }' N; \  v& q  r- G5 J

可以发现只要通过参数设定就可以将所要求解的方程转化成这种标准形式。

具体求解步骤如下：

(1)  设置参数

• 
  
  c = 1;
  , a: d( ?! M( f. P
  
• 
  1 U7 Z- V7 X& @( B4 o5 d; H, q
  6 y: D5 u7 T' B
  # P% s8 {% n1 W; E7 X6 n% s
  
• 
  % @* J/ }- l6 m1 x( h
  a = 0;
  $ F/ I: W5 H! l6 o5 u5 a8 O3 y
  ) H+ r# z( y: l$ P8 B4 ~7 y" ]
• 
  $ }& V- i2 a2 W" M
  ( L8 u9 ?: _: l. m; N( j
  
  
• 
  
  ! O2 ^% I* e- Gf = 0;
  2 L9 j1 c% Y5 Z' F; a
  
• 
  , g4 e7 d* F# Q4 w8 ]5 g
  
  ) G/ z  u7 @% j
  * [7 W* f$ `$ i
• 
  
  & y( V9 {* E2 A, nm = 1;
  
  5 `  d  A! I& m3 m. E9 Q6 ~5 m
  

(2)  定义波的空间位置

• 
  
  2 Z/ ?" b8 F9 MnumberOfPDE = 1;
  
  2 v$ X) E: y+ F+ v
• 
  8 b* G3 o5 u; W: r: @, g, `
  , _2 }8 C) ~7 Z9 Y9 g0 Z
  : }/ Y7 L  w- t
  / K7 p$ W: `& n3 r- i1 C
• 
  + K0 p: F& A/ }) W
  $ E/ D" |' l0 V' h5 B+ t8 ~model = createpde(numberOfPDE);
  ; t) t8 k) ]* x7 |& H9 e
  ( e+ t* k3 R5 Z' P8 P9 |& }
• 
  
  
  * v5 l6 Z  ]- r; ^9 ^
  
• 
  
  0 J6 m, v  e6 Q7 Z6 n$ OgeometryFromEdges(model,@squareg);
  & A: P0 r+ a  \. D
  
• 
  
  
  
  
• 
  
  ) ~# j/ s% J" R$ i5 P% h! E4 Epdegplot(model,'EdgeLabels','on');
  5 X6 W& T7 j: R: W& |9 w
  7 r9 a( Y+ m) i" f0 J0 T
• 
  & [1 x5 t4 W) f4 O- A! A
  " l8 O8 I7 |2 Z4 C' g* K; |
  
  0 \( `$ o0 t1 I/ M. i
• 
  % e+ I0 N  x$ r$ G' K7 C( q
  ylim([-1.1 1.1]);
  ) C+ w) s# J9 I
  ( G. |; k% _5 B
• 
  6 i1 l$ w4 z+ q7 v
  
  + |3 `) P* g) {. u* O3 L4 R" ?
  
• 
    H. x' ^6 n) ~5 v: m+ _- S7 C, f
  axis equal
  
  
• 
  
  
  1 a$ ~( J' g" s; T( C( x
  % [8 f& S8 O8 e% w/ i
• 
  8 p2 J( e4 @4 D" h. M0 b4 o
  ' z% x# h/ w' K1 ]4 _title 'Geometry With Edge Labels Displayed';
  ) ^! i5 {: I6 c, I( p
  
• 
  , x$ i9 A& r/ V% n- Y. z3 V& C
  ; t4 p1 U1 W& |% ~3 p% x+ ]3 w
  0 K. }) b4 z8 p. Q
  6 D8 w  I6 ^2 {8 U5 Y1 a0 _; B
• 
  ' k1 P. Q" R) i4 _
  . g5 |4 C5 }* f0 i* r+ Hxlabel x
  % L: A0 X/ C8 @% A/ m0 ^
  ; s) ^( U1 o' d" x
• 
  % z; v3 m- R9 A4 _# L( ]+ u
  
  
  
• 
  
  ) k6 M! j* Z3 R* W& ~, Xylabel y
  
  
  

(3)  定义微分方程模型的系数和边界条件

• 
  9 W/ X! `2 ^' j; z( \
  6 c7 P1 ?) T) x# I' DspecifyCoefficients(model,'m',m,'d',0,'c',c,'a',a,'f',f);
  
  
• 
  
  
  9 F, U6 l4 ^6 A; _
  
• 
  
  - o# q+ f: t* c2 b7 z. tapplyBoundaryCondition(model,'dirichlet','Edge',[2,4],'u',0);
  1 k! Q3 t/ C# g9 d* J/ X; Y
  
• 
  0 S& S( e& h9 e% l' j
  
  
  
• 
  
  applyBoundaryCondition(model,'neumann','Edge',([1 3]),'g',0);
  
  
  

4)  定义该问题的有限元网格

• 
  6 ?" l/ Q) f3 @7 m
  8 {, t9 F( Z1 Y0 V2 xgenerateMesh(model);
  
  
• 
  
  3 G3 B0 i* x4 {
  
  " D- E0 v. b, D
• 
  * M+ @7 R: h% `
  figure
  
  
• 
  
  8 `' g" z0 r( U0 v
  0 @- }& n' R5 V! ^
  ; i) b# f, p" I0 ?1 U0 _5 z( p5 w/ Q2 \0 J
• 
  ' L8 l' x2 c2 y% M8 P
  pdemesh(model);
  
  
• 
  - Z# F% P, y. F+ B9 S
  
  
  / h/ x2 N2 U0 \4 ]- q6 o( ~
• 
  
  ylim([-1.1 1.1]);
  . V1 n6 P: S5 a" c6 {* g2 B
  
• 
  7 w; H# p2 s* g) m' d2 _6 A
  0 L* ?, S7 N" K% k9 w9 e# X, ~
  
  
• 
  9 ]/ R4 a/ K5 Y3 o" N
  axis equal
  
  8 T5 I; m( S+ z0 G! t
• 
  0 K3 Y! o- O, I0 B
  4 G4 r# \2 P7 x; y1 O
  
  
• 
  % F1 H: H1 h4 u8 I! {" y4 a
  xlabel x
  0 @. Y6 i+ \) e: [, @3 o
  
• 
  
  1 H0 k& P1 q  ]# k' \1 s
  
  
• 
  
  % `; R+ q1 x1 Z4 H7 M! C& dylabel y
  
  % `# I2 e" r: D, ]7 h" Q
  : V; [& {3 ]! C$ v# P+ ~3 Q

(5)  定义初始条件

• 
  
  u0 = @(location) atan(cos(pi/2*location.x));
  
  $ g! a5 h, l; b! a+ ]1 O4 @$ c) o5 k
• 
  8 v& N' P( g" m7 r
  ( n+ W/ m" ^2 P- ^8 p
  
  7 Z& g) m$ A0 U/ V! q6 V/ t7 S
• 
  
  ut0 = @(location) 3*sin(pi*location.x).*exp(sin(pi/2*location.y));
  
  
• 
  
  ) I8 `- v  `. K% _% @$ P& I
  
  % C4 [6 U& ~: z4 t7 d( y
• 
  # ]: `. ^: J: H2 A0 Q# I
  ) E7 O: V% _* @* c" `setInitialConditions(model,u0,ut0);
  
  # Y9 _$ E$ ^. D3 i" ?7 J
  

(6)  方程的求解

• 
  $ q/ `, c; G0 [% i6 T. G
  & N, Z. L4 G7 Yn = 31;  % 求解次数
  
  " |+ `7 K! [; \+ s
• 
  
  & O% L) Y8 F1 |# Q; m/ ~9 [
  & X- W9 w0 \. T0 g! ^' ]
  ) K0 z6 h3 |) M+ J' J6 R
• 
  
  ' O4 L- s& [3 Y' Ntlist = linspace(0,5,n);
  
  6 ^( M7 s( [5 O3 ~1 r, @
• 
  ) [+ p1 Z' ]5 R/ _* B2 V. ~1 u3 u) F
  2 Q1 C* V2 ?0 K9 }% ]. S* Q
  
  
• 
  1 ^- m0 l- M# X4 p  u
  / H& ]! K5 }5 ?  t# W- vmodel.SolverOptions.ReportStatistics ='on';
  . O$ z7 @( E6 F6 V) d8 N7 @/ o
  : Q9 n2 [. r# A0 M4 o5 p
• 
  
  
  
  
• 
  
  result = solvepde(model,tlist);
  
  
• 
  
  
  
  
• 
  % X, W3 c: ~3 j/ h. a; c
  - i9 y6 c4 v( m4 vu = result.NodalSolution;
  
  ) }$ F  d' D( E0 x9 o
    b4 y6 ?5 X7 s+ q$ u' ?% m

(7)  模型的数值仿真

• 
  * K3 S- Q* O3 D0 j  {8 |) E5 x
  figure
  
  
• 
  0 ]& `. `* ]" R# i1 P$ H2 G
  ) i+ z* y4 Y: s3 |7 t
  / J( s$ Z  X7 a+ Q& a( u  u' V
  
• 
  
  umax = max(max(u));
  
  
• 
  
  
  ' w8 M' m- L" ^( [
  
• 
  
  ; ]9 M& c7 q: m, G, C# C# R5 m5 @umin = min(min(u));
  ' H( {; w; {( Z) T
  
• 
  5 J) I! n1 c2 E/ {' t
  : f5 }( h8 Y& G, Z. N8 d
  5 ~8 M( e# D0 P7 L: E' Y, ?( j
  8 y! i& ?% c5 b3 X3 ~6 k4 r
• 
  
  for i = 1:n
  
  
• 
  ! z. y( z3 D  ]  o  v( \  |1 F
  * a( y# l" T5 C
  - \* J5 C/ p- s: E+ ]1 N: T8 j
  $ |2 {' h8 `- B
• 
  
      pdeplot(model,'XYData',u(:,i),'ZData',u(:,i),'ZStyle','continuous',...
  
  
• 
  
  
  ( h7 V! ?5 n% b9 K' v
  
• 
  8 T7 N4 ?+ }$ c" }
                    'Mesh','off','XYGrid','on','ColorBar','off');
  
  % R( N* i& v3 ^5 ^, ~; k
• 
  ) f# L8 B( C% O8 m! _- a
  ; ]3 k5 r( ?& m6 [
  ; d" N  o! {. u+ N1 Q3 Y  h" U
  
• 
  2 F8 l8 j9 d0 Y) D8 `$ b: e% @: B
      axis([-1 1 -1 1 umin umax]);
  
  8 j* R* d4 M/ C( ~( o5 j3 K2 P
• 
  0 s; }4 U8 c9 z9 v" b, Y: y$ C
  , m. J$ j5 V& K7 a& }+ @$ h
  + x2 N  Q3 V; r2 v5 b# w
  6 j+ Q" Y$ `0 D
• 
  
  6 U3 _8 Y* {  _' Z    caxis([umin umax]);
  
  ' }+ F* u, n: R+ }
• 
  
  
  
  
• 
  ) o+ n* o# @0 p$ r: s
      xlabel x
  0 F- x2 J! A" n- b+ b# l: {
  
• 
  
  % n* w. a: s+ T/ a1 ]. s4 u
  6 B9 R4 T  ^3 i5 q
  
• 
  : v" [9 [7 @( k/ p/ T- r
  8 z6 W) z" `/ g6 b2 |    ylabel y
  " v* Z: }. B# p# U% ~) g
  
• 
  
  2 m( a! b4 |) F$ j$ P4 m
  
  5 h3 H+ [! ?4 j, z: w# l
• 
  
  3 m' O7 o. Z8 K- T+ Q1 [* X    zlabel u
  ' U6 j3 [4 k: [) v) j" s# V
  
• 
  9 m0 x  b$ n$ k% E
  7 O6 V7 P! q0 m8 c& c2 w
  
  
• 
  % _/ I3 {6 q$ l! S, Q
      M(i) = getframe;
  
  
• 
  
  # V) z2 A& r% L. o
  ; V8 f+ Z  y5 V
  - x! [: p- D* l4 D
• 
  
  ; x# N/ S8 m9 J8 Lend
  
  ! u6 R8 w/ \$ o& [0 C+ R2 _
• 
  
  0 _2 x8 R& J1 g+ p5 j) D/ C. S
  
  6 y* j, t+ {: v: `5 D6 }+ K# L
  / ~! M/ R5 n3 V( j" h6 s2 ?

9 |) ~! T# ]4 A- V: _, q  E& K
! Q( R, b3 B' `8 F