例如，我们求下列方程组在
(
1
,
1
,
1
)
(1,1,1)
(1,1,1) 附近的解并对结果进行验证。
{
sin
x
+
y
+
z
2
e
x
=
0
x
+
y
+
z
=
0
x
y
z
=
0
\left\{\begin{matrix}\sin x+y+z^{2}e^{x}=0 \\x+y+z=0 \\xyz=0 \end{matrix}\right.
⎩⎨⎧sinx+y+z2ex=0x+y+z=0xyz=0
首先，我们建立函数文件 myfun.m。
function F =myfun(X)
x=X(1);
y=X(2);
z=X(3);F(1)=sin(x)+y+z^2*exp(x);F(2)=x+y+z;F(3)=x*y*z;
end
在给定的初值
x
0
=
1
，
y
0
=
1
，
z
0
=
1
x_{0}=1，y_{0}=1，z_{0}=1
x0=1，y0=1，z0=1 下，调用 fsolve 函数求方程的根，程序如下：
>> option=optimset('Display','off');>> X=fsolve(@myfun,[1,1,1],option)
X =0.0224-0.0224-0.0000
将求得的解代回原方程，可以验证结果是否正确，程序如下：
>> q=myfun(X)
q =1.0e-06*-0.5931-0.00000.0006
可见得到了较高精度的结果。
例如，我们求圆和直线的两个交点。
圆：
x
2
+
y
2
+
z
2
=
9
圆：x^{2}+y^{2}+z^{2}=9
圆：x2+y2+z2=9
直线：
{
3
x
+
5
y
+
6
z
=
0
x
−
3
y
−
6
z
−
1
=
0
直线：\left\{\begin{matrix}3x+5y+6z=0 \\x-3y-6z-1=0 \end{matrix}\right.
直线：{3x+5y+6z=0x−3y−6z−1=0
该问题即为求解方程组
{
x
2
+
y
2
+
z
2
=
9
3
x
+
5
y
+
6
z
=
0
x
−
3
y
−
6
z
−
1
=
0
\left\{\begin{matrix}x^{2}+y^{2}+z^{2}=9 \\3x+5y+6z=0 \\x-3y-6z-1=0 \end{matrix}\right.
⎩⎨⎧x2+y2+z2=93x+5y+6z=0x−3y−6z−1=0
无约束最优化问题的一般描述为
m
i
n
x
f
(
x
)
\underset{x}{min} f(x)
xminf(x)
其中，
x
=
[
x
1
,
x
2
,
…
,
x
n
]
T
x=[x_{1},x_{2},…,x_{n}]^{T}
x=[x1,x2,…,xn]T，该数学表示的含义亦即求取一组
x
x
x，使得目标函数
f
(
x
)
f(x)
f(x) 为最小，时间最短等。MATLAB 提供了 3 个求最小值的函数，它们的调用格式如下。
（1） [x,fval]=fminbnd(filename,x1,x2,option)：求一元函数在
(
x
1
,
x
2
)
(x1, x2)
(x1,x2) 区间中的极小值点
x
x
x 和最小值 fval。
（2） [x,fval]=fminsearch(filename,x0,option): 基于单纯形算法求多元函数的极小值点
x
x
x 和最小值 fval。
（3） [x,fval]=fminunc(ilename,x0,option)：基于拟牛顿法求多元函数的极小值点
x
x
x 和最小值 fval。
MATLAB 没有专门提供求函数最大值的函数，但可以注意到
−
f
(
x
)
-f(x)
−f(x) 在区间
(
a
,
b
)
(a, b)
(a,b) 上的最小值就是
f
(
x
)
f(x)
f(x) 在
(
a
,
b
)
(a, b)
(a,b) 的最大值， 所以 fminbnd(- f,x1,x2) 返回函数
f
(
x
)
f(x)
f(x) 在区间
(
x
1
,
x
2
)
(x1, x2)
(x1,x2) 上的最大值。
例如，我们求
f
(
x
)
=
x
−
1
x
+
5
f(x)=x-\frac{1}{x}+5
f(x)=x−x1+5 在区间
(
−
10
,
−
1
)
(-10, -1)
(−10,−1) 和
(
1
,
10
)
(1, 10)
(1,10) 上的最小值点。
程序如下：
>> f=@(x) x-1./x+5;>>[x,fmin]=fminbnd(f,-10,-1)%求函数在(-10，-1)内的最小值点和最小值
x =-9.9999
fmin =-4.8999>>fminbnd(f,1,10)%求函数在(1,10)内的最小值点
ans =1.0001
例如，设
f
(
x
,
y
,
z
)
=
x
+
y
2
4
x
+
2
z
+
z
2
y
f(x,y,z)=x+\frac{y^{2}}{4x}+\frac{2}{z}+\frac{z^{2}}{y}
f(x,y,z)=x+4xy2+z2+yz2 我们求函数
f
f
f 在
(
0.5
,
0.5
,
0.5
)
(0.5, 0.5, 0.5)
(0.5,0.5,0.5) 附近的最小值点。
首先，我们建立函数文件 fxyz0.m。
function f=fxyz0(u)
x=u(1);
y=u(2);
z=u(3);
f=x+y.^2./x/4+z.^2./y+2./z;
end
随后，在 MATLAB 命令行窗口输入以下程序：
>>[X,fmin]=fminsearch(@fxyz0,[0.5,0.5,0.5])
X =0.50001.00001.0000
fmin =4.0000
2. 有约束最优化问题求解
有约束最优化问题的一般描述为
m
i
n
x
s
.
t
.
G
(
x
)
≤
0
f
(
x
)
\underset{x s.t. G(x)≤0}{min} f(x)
xs.t.G(x)≤0minf(x)
其中，
x
=
[
x
1
,
x
2
,
…
,
x
n
]
T
x=[x_{1},x_{2},…,x_{n}]^{T}
x=[x1,x2,…,xn]T，该数学表示的含义亦即求取一组
x
x
x，使得目标函数
f
(
x
)
f(x)
f(x) 为最小，且满足约束条件
G
(
x
)
≤
0
G(x)≤0
G(x)≤0。记号 s.t. 是英文 subject to 的缩写，表示
x
x
x 要满足后面的约束条件。
约束条件可以进一步细化如下。
（1） 线性不等式约束：
A
x
≤
b
Ax≤b
Ax≤b。
（2） 线性等式约束：
A
e
q
x
=
b
e
q
A_{eq}x=b_{eq}
Aeqx=beq。
（3） 非线性不等式约束：
C
x
≤
0
Cx≤0
Cx≤0。
（4） 非线性等式约束：
C
e
q
x
=
0
C_{eq}x=0
Ceqx=0。
（5）
x
x
x 的下界和上界：
L
b
n
d
≤
x
≤
U
b
n
d
L_{bnd}≤x≤U_{bnd}
Lbnd≤x≤Ubnd。
例如，我们求解有约束最优化问题。
m
i
n
x
s
.
t
.
{
x
1
+
0.5
x
2
≥
0.4
0.5
x
1
+
x
2
g
e
0.5
x
1
≥
0
,
x
2
≥
0
f
(
x
)
=
0.4
x
2
+
x
1
2
+
x
2
2
−
x
1
x
2
+
1
30
x
1
3
\underset{x s.t.\left\{\begin{matrix}x_{1}+0.5x_{2}\ge 0.4 \\0.5x_{1}+x_{2}ge 0.5 \\x_{1}\ge 0,x_{2}\ge 0 \end{matrix}\right.}{min}f(x)=0.4x_{2}+x_{1}^{2}+x_{2}^{2}-x_{1}x_{2}+\frac{1}{30}x_{1}^{3}
xs.t.{x1+0.5x2≥0.40.5x1+x2ge0.5x1≥0,x2≥0minf(x)=0.4x2+x12+x22−x1x2+301x13
首先，我们编写定义目标函数的函数文件 fop.m。
function f=fop(x)
f=0.4*x(2)+x(1)^2+x(2)^2-x(1)*x(2)+1/30*x(1)^3;
end
线性规划是研究线性约束条件下线性目标函数的极值问题的数学理论和方法。线性规划问题的标准形式为
m
i
n
x
s
.
t
.
{
A
x
=
b
A
e
q
x
=
b
e
q
L
b
n
d
≤
x
≤
U
b
n
d
f
(
x
)
\underset{x s.t.\left\{\begin{matrix}Ax=b \\A_{eq}x=b_{eq} \\L_{bnd}≤x≤U_{bnd} \end{matrix}\right.}{min} f(x)
xs.t.{Ax=bAeqx=beqLbnd≤x≤Ubndminf(x)
例如，我们求解线性规划问题。
m
i
n
x
s
.
t
.
{
x
1
−
x
2
+
x
3
≤
20
3
x
1
+
2
x
2
+
4
x
3
≤
42
3
x
1
+
2
x
2
≤
30
x
1
≥
0
,
x
2
≥
0
,
x
3
f
(
x
)
\underset{x s.t.\left\{\begin{matrix}x_{1}-x_{2}+x_{3}\le 20 \\3x_{1}+2x_{2}+4x_{3}\le 42 \\3x_{1}+2x_{2}\le 30 \\x_{1}\ge 0,x_{2}\ge 0,x_{3} \end{matrix}\right.}{min} f(x)
xs.t.⎩⎨⎧x1−x2+x3≤203x1+2x2+4x3≤423x1+2x2≤30x1≥0,x2≥0,x3minf(x)
考虑常微分方程的初值问题：
y
′
=
f
(
t
,
y
)
,
t
0
≤
t
≤
T
y'=f(t,y),t_{0}\le t\le T
y′=f(t,y),t0≤t≤T
y
(
t
0
)
=
y
0
y(t_{0})=y_{0}
y(t0)=y0
所谓其数值解法，就是求它的解
y
(
t
)
y(t)
y(t) 在节点
t
0
<
t
1
<
⋯
<
t
m
t_{0}<t_{1}<\dots < t_{m}
t0<t1<⋯<tm 处的近似值
y
0
,
y
1
,
…
,
y
m
y_{0},y_{1},\dots ,y_{m}
y0,y1,…,ym 的方法。所求得的
y
0
,
y
1
,
…
,
y
m
y_{0},y_{1},\dots ,y_{m}
y0,y1,…,ym 称为常微分方程初值问题的数值解。
一般采用等距节点
t
n
=
t
0
+
n
h
,
n
=
0
,
1
,
…
,
m
t_{n}=t_{0}+nh, n=0, 1, …,m
tn=t0+nh,n=0,1,…,m，其中
h
h
h 为相邻两个节点间的距离，叫做步长。
对于一阶常微分方程的处置问题，在求解未知函数
y
y
y 时，
y
y
y 在
t
0
t_{0}
t0 点的值
y
(
t
0
)
=
y
0
y(t_{0})=y_{0}
y(t0)=y0 是已知的，并根据高等数学中的中值定理，应有
{
y
(
t
0
+
h
)
=
y
1
≈
y
0
+
h
f
(
t
0
,
y
0
)
y
(
t
0
+
2
h
)
=
y
2
≈
y
1
+
h
f
(
t
1
,
y
1
)
,
h
>
0
\left\{\begin{matrix}y(t_{0}+h)=y_{1}\approx y_{0}+hf(t_{0},y_{0}) \\y(t_{0}+2h)=y_{2}\approx y_{1}+hf(t_{1},y_{1}) \end{matrix}\right.,h>0
{y(t0+h)=y1≈y0+hf(t0,y0)y(t0+2h)=y2≈y1+hf(t1,y1),h>0
一般地，在任意点
t
i
=
t
0
+
i
h
t_{i}=t_{0}+ih
ti=t0+ih，有
y
(
t
0
+
i
h
)
=
y
i
≈
y
i
−
1
+
h
f
(
t
i
−
1
,
y
i
−
1
)
,
i
=
1
,
2
,
…
,
n
y(t_{0}+ih)=y_{i}\approx y_{i-1}+hf(t_{i-1},y_{i-1}),i=1,2,…,n
y(t0+ih)=yi≈yi−1+hf(ti−1,yi−1),i=1,2,…,n
当
(
t
0
,
y
0
)
(t_{0},y_{0})
(t0,y0) 确定后，根据上述递推式能计算出未知函数
y
y
y 在点
t
i
=
t
0
+
i
h
,
i
=
1
,
…
,
n
t_{i}=t_{0}+ih,i=1,…,n
ti=t0+ih,i=1,…,n 的一列数值解：
y
i
=
y
1
,
y
2
,
…
,
y
n
,
i
=
1
,
2
,
…
,
n
y_{i}=y_{1},y_{2},…,y_{n},i=1,2,…,n
yi=y1,y2,…,yn,i=1,2,…,n
当然，地推过程中有一个误差累积的问题。在实际计算过程中，使用的递推公式一般进行过改造，著名的龙格——库塔公式是
y
(
t
0
+
i
h
)
=
y
i
≈
y
i
−
1
+
h
6
(
k
1
+
2
k
2
+
2
k
3
+
k
4
)
y(t_{0}+ih)=y_{i}\approx y_{i-1}+\frac{h}{6}(k_{1}+2k_{2}+2k_{3}+k_{4})
y(t0+ih)=yi≈yi−1+6h(k1+2k2+2k3+k4)
其中
k
1
=
f
(
t
i
−
1
,
y
i
−
1
)
k_{1}=f(t_{i-1},y_{i-1})
k1=f(ti−1,yi−1)
k
2
=
f
(
t
i
−
1
+
h
2
,
y
i
−
1
+
h
2
k
1
)
k_{2}=f(t_{i-1}+\frac{h}{2},y_{i-1}+\frac{h}{2}k_{1})
k2=f(ti−1+2h,yi−1+2hk1)
k
3
=
f
(
t
i
−
1
+
h
2
,
y
i
−
1
+
h
2
k
2
)
k_{3}=f(t_{i-1}+\frac{h}{2},y_{i-1}+\frac{h}{2}k_{2})
k3=f(ti−1+2h,yi−1+2hk2)
k
4
=
f
(
t
i
−
1
+
h
,
y
i
−
1
+
h
k
3
)
k_{4}=f(t_{i-1}+h,y_{i-1}+hk_{3})
k4=f(ti−1+h,yi−1+hk3)
2. 龙格—库塔法的实现
MATLAB 提供了多个求常微分方程初值问题数值解的函数，一般调用格式如下：
[t,y]=solver(filename,tspan,y0)
其中，
t
t
t 和
y
y
y 分别给出时间向量和相应的状态向量。solver 为求常微分方程数值解的函数，下表中列出了常用函数采用的方法和适用的场合。filename 是定义
f
(
t
,
y
)
f(t, y)
f(t,y) 的函数，该函数必须返回一个列向量。tspan 形式为
[
t
0
,
t
f
]
[t0, tf]
[t0,tf]，表示求解区间。
y
0
y0
y0 是初始状态列向量。
例如，已知一个二阶线性系统的微分方程为
{
d
2
x
d
t
2
+
a
x
=
0
,
a
>
0
x
(
0
)
=
0
,
x
′
(
0
)
=
1
\left\{\begin{matrix}\frac{\mathrm{d}^{2} x}{\mathrm{d} t^{2}}+ax=0,a>0 \\x(0)=0,x'(0)=1 \end{matrix}\right.
{dt2d2x+ax=0,a>0x(0)=0,x′(0)=1 取
a
=
2
a=2
a=2，绘制系统的时间响应曲线和相平面图。
函数 ode23 和 ode45 是对一阶常微分方程组设计的，因此对高阶常微分方程，需先将它转化为一阶常微分方程组，即状态方程。令
x
2
=
x
，
x
1
=
x
′
x_{2}=x，x_{1}=x'
x2=x，x1=x′， 则得到系统的状态方程：
{
x
2
′
=
x
1
x
1
′
=
−
a
x
2
x
2
(
0
)
=
0
,
x
1
(
0
)
=
1
\left\{\begin{matrix}x'_{2}=x1 \\x'_{1}=-ax_{2} \\x_{2}(0)=0,x_{1}(0)=1 \end{matrix}\right.
⎩⎨⎧x2′=x1x1′=−ax2x2(0)=0,x1(0)=1
首先，我们建立函数文件 sys.m。
function xdot=sys(t,x)
xdot=[-2*x(2);x(1)];
end
取
t
0
=
0
，
t
f
=
20
t0=0，tf=20
t0=0，tf=20，求微分方程的解，程序如下：
>>[t,x]
ans =01.000000.00011.00000.00010.00011.00000.00010.00021.00000.00020.00021.00000.0002
…
19.6681-0.89680.311619.7511-0.94220.235219.8340-0.97470.155619.9170-0.99370.073820.0000-0.9991-0.0090
结果第一列为
t
t
t 的采样点，第二列和第三列分别为
x
′
x'
x′ 和
x
x
x 与
t
t
t 对应点的值。
