浅谈 MATLAB R2014b solve 函数新功能

winner245

最近发现 2014b solve 函数功能有所加强，比如：一些之前解不了的方程已经可以求解了，一些求不到符号解的方程会自动转向数值求解。
不过，比较显著的变化是：1. 对于解集含参变量的支持更好了， 2. 对于不等式（组）的求解有更好的支持了

1. 对于解集含参变量的支持更好了

这是通过新增的 ReturnConditions 参数调用实现的：

S = solve(eqns,vars,'ReturnConditions',true)

这一功能的引入，使得 solve 函数可以针对解集中含有参变量的情形做更好的处理，返回的 S 将是一个结构体 （当然也可以指定多个独立的返回变量），由三部分组成，分别是：x  (求解变量)、parameters （参变量）、conditions（参变量满足的条件，即解存在的条件）。唯一需要做的是将这里的 ReturnConditions 设置为 true （默认是 false）。这一功能的引入，使得求解方程（组）的通解表达式变得方便很多。

这里以方程中含周期函数为例。一般，周期函数方程的解集需要以通解的形式给出，通解表达式中的参变量所满足的条件可以由ReturnConditions参数来指定返回。而在之前的版本里，要么干脆得不到通解而只能得到若干特解，要么得到一种十分不友好的通解表达式。下面分类示例：

首先 sin(x) = 1/2 为例，如果我们直接写下面的代码：

>> syms x
>> solve(tan(x)-1/2)

ans =
atan(1/2)

则无论是2014b 还是之前的版本都只能得到以上一个特解。在 2014b 之前，我们似乎除了上述求解方法，也没有其他快捷的方法 “直接求得” 通解表达式了。但 2014b 中却可以通过ReturnConditions 轻易得到通解了：

>> syms x
>> [x, p, c] = solve(tan(x)-1/2,'ReturnConditions',true)

x =
atan(1/2) + pi*k

p =
k

c =
in(k, 'integer')

可以看到，我们得到了通解表达式，而且通解表达式中参变量 k 满足的条件也一并返回到了。

再举一个例子：

>> syms x
>> solve(sin(x)^2 + 2*sin(x) - 1)

ans =
     -asin(1 - 2^(1/2))
     -asin(2^(1/2) + 1)
pi + asin(1 - 2^(1/2))
pi + asin(2^(1/2) + 1)

2014b 中得到的是以上4个特解。这段程序如果放在老一点的版本里，比如 2014a 或 2012b 里，会得到下列的结果：
ans =
Z_ intersect (Dom::ImageSet(pi + asin(2^(1/2) + 1) + 2*pi*k, k, Z_) union Dom::ImageSet(- asin(2^(1/2) + 1) + 2*pi*k, k, Z_) union Dom::ImageSet(pi + asin(- 2^(1/2) + 1) + 2*pi*k, k, Z_) union Dom::ImageSet(- asin(- 2^(1/2) + 1) + 2*pi*k, k, Z_))

这实际上是一个通解表达式，可是这个表达式多少看起来让人不舒服，恐怕一般人很不习惯读这样的表达式吧，而且，即使用了pretty指令，也得不到改观，依然很难读。那么为了得到更加 user-friendly 的solution，我们看看 2014b 给我们带来什么样的改变：

>> syms x
>> [x, p, c] = solve(sin(x)^2 + 2*sin(x) - 1,'ReturnConditions',true)

x =
      2*pi*k - asin(1 - 2^(1/2))
      2*pi*k - asin(2^(1/2) + 1)
pi + asin(1 - 2^(1/2)) + 2*pi*k
pi + asin(2^(1/2) + 1) + 2*pi*k

p =
k

c =
in(k, 'integer')
in(k, 'integer')
in(k, 'integer')
in(k, 'integer')

可以看到4个通解表达式中参变量 k 所满足的条件都一一给出了，而且通解表达式非常友好！

其实，这里我主要想说的是，ReturnConditions 最大的好处是，它把通解集中参变量满足的条件以表达式的形式返回给conditions参数了，利用 conditions 返回值，我们还可以做一系列后续计算。比如，我们根据这个 conditions 返回值，只需结合 subs、isAlaways 函数做一些极其简单的操作，便可以对特定的参变量值做判断，从而筛选出我们感兴趣的某些解，并且后续一系列计算变得很方便。而早期的版本的条件表达式并没有直接返回给我们，要做上述操作，我们还需要去手写这些表达式，这是极其不便利的。下面依然回到第一个例子说明：

1. syms x
   
2. [x, p, c] = solve(tan(x)-1/2,'ReturnConditions',true);
   
3. k = 15;
   
4. if isAlways(subs(c, p, k))
   
5.     x = subs(x,p,k)
   
6. end

复制代码

x =
15*pi + atan(1/2)

显然，以上操作正是基于返回得到的条件表达式做了后续计算，而在早起版本里，这些都不是直接能做到的。

2. 对于求解不等式（组）的支持也更好了

例子 1：求 [-2*pi, 2*pi] 之间所有自然对数小于1.5的整数

>> syms x
>> assume(x,'integer')
>> solve(x<=2*pi,x>=-2*pi,log(x)<1.5)

ans =
  1
  2
  3
  4

而相同的程序在之前的版本中都会给出 Warning: 3 equations in 1 variables. 的警告，因为它把3个不等式看作3个方程了，而这里只有一个未知数 x。

需要说明的是，上面为了力求说明问题，例子都很简单。下面，有时间我再继续补充

winner245

发现这个增强版的 solve 函数也有不尽如人意的地方。比如，求解方程：x+x^0.4-5=0。先看看三种调用形式结果的差别

1. 最原始的 string 方程形式 （注意，MATHWORKS 已经明确宣布 string 形式的方程在今后的版本中将不再推荐或弃用 string 方程，因为其局限性十分明显）

1. solve('x+x^.4-5')

复制代码

ans =

3.3735683862761403832267961226935

2. 较常见的符号（sym）方程输入方式

1. syms x
   
2. solve(x+x^.4-5)

复制代码

Warning: The solutions are parameterized by the symbols: z2. To include parameters and conditions in the solution,
specify the 'ReturnConditions' option.
> In solve>warnIfParams at 511
  In solve at 361
  In Untitled at 15
Warning: The solutions are valid under the following conditions: z2^2 + z2^5 - 5 == 0 & angle(z2) in
Dom::Interval(-pi/5, [pi/5]) == 0. To include parameters and conditions in the solution, specify the
'ReturnConditions' option.
> In solve>warnIfParams at 518
  In solve at 361
  In Untitled at 15

ans =

z2^5

得到了参数解。而且推荐我们使用 ReturnConditions 选项，于是我们使用ReturnConditions来求解如下

3. 使用ReturnConditions 方式，根绝返回条件筛选出解

1. syms x
   
2. S = solve(x+x^.4-5,'ReturnConditions',true);
   
3. S.x
   
4. S.parameters
   
5. S.conditions
   
6. syms z2
   
7. z3 = solve(z2^5 + z2^2 - 5);
   
8. id = isAlways(subs(S.conditions, S.parameters, z3));
   
9. vpa(subs(S.x,S.parameters,z3(id)))

复制代码

ans =
z2^5

ans =
z2^5

ans =
z2

ans =
z2^5 + z2^2 - 5 == 0 & angle(z2) in Dom::Interval(-pi/5, [pi/5]) == 0

ans =
3.3735683862761403832267961226935

虽然最终成功筛选到了解，可是这种方式明显比第 1 种方法 solve('x+x^.4-5') 复杂了好多，这是不太理解的地方。因为很多复杂的场合，我们显然是一开始 syms 到底，不会使用 string 形式的方程，况且这也不是mathworks推荐的方式



注： 第1种方法solve('x+x^.4-5')在之前的版本（比如 2012b）里会产生错误解：

>> x = solve('x+x^.4-5=0')
>> x+x.^.4-5

x =
                                       3.3735683862761403832267961226935
6.6276037792645952670570079999919 + 1.4041584183916010600989318333578*i
6.6276037792645952670570079999919 - 1.4041584183916010600989318333578*i
4.1856120275973345413295939386613 - 1.6319622743772952871024038130991*i
4.1856120275973345413295939386613 + 1.6319622743772952871024038130991*i

ans =
                                    1.9601368299961644195378681758181e-36
  3.7697065070114507754018039893727 + 1.5834649436378311725014255197239*i
  3.7697065070114507754018039893727 - 1.5834649436378311725014255197239*i
0.98936010647869987152438237429449 - 1.9021871853618729040016942529236*i
0.98936010647869987152438237429449 + 1.9021871853618729040016942529236*i

rocwoods

winner245 发表于 2014-11-6 01:30
发现这个增强版的 solve 函数也有不尽如人意的地方。比如，求解方程：x+x^0.4-5=0。先看看三种调用形式结果 ...

感谢刘兄的测试工作！MATLAB的符号计算的确相对另外两个数学软件Maple和MMTC偏弱，在MATLAB里，个人建议符号计算还是作为辅助比较好，我一般都是用它推一推公式，然后用matlabFunction转成句柄再数值求解。

winner245

rocwoods 发表于 2014-11-6 09:20
感谢刘兄的测试工作！MATLAB的符号计算的确相对另外两个数学软件Maple和MMTC偏弱，在MATLAB里，个人建议 ...

的确，在符号计算领域 MATLAB 还是逊色很多。不过尤其是最近几次MATLAB更新，我都对它寄予了很高的期待，以至于会迫不及待的感受一下新版的一些改进功能

winner245

2014b solve 求解不等式变得很不方便

先看看极其简单的不等式例子：
>> syms x
>> solve(x^2-1>0)

各个版本结果如下：

2012b 和 2014a 得到

ans =
   Dom::Interval(1, Inf)
Dom::Interval(-Inf, -1)

而 2014b 却给出：
ans =
-2
  2
这个结果让人有点莫不着头脑。

再比如一个最简单的不等式：

>> syms x
>> solve(x-3<0)

2014b 居然给出的结果是：
ans =
2

而 2014a 和 2012b 的结果是
ans =
Dom::Interval(-Inf, 3)

也不知道 2014b 怎么回事，默认的 solve 调用方式根本得不到不等式的解（也不知道这是不是 bug）。不过，2014b 可以通过设置 ReturnConditions 为true来实现求解：

>> syms x
>> S = solve(x^2-1>0,x,'ReturnConditions',true);
>> S.conditions

ans =

x < -1
  1 < x

rocwoods

winner245 发表于 2014-11-13 11:53
2014b solve 求解不等式变得很不方便

先看看极其简单的不等式例子：

以下只是个人猜测，默认情况下，2014b解不等式是给出解区间的一个特解，有几个区间就给几个特解。如果要返回全部解就按照刘兄说的那样设置。
至于为什么要这么改，我就想不明白了。