求推荐可代替 matlab 部分功能的第三方数学库

[kindlebeta]

小弟最近在编一个工程方面的计算软件，在求解非线性方程组的根时遇到了一些问题，特来向各位大神请教。
方程组如下所示（双曲线方程）， s,b,h2 未知量，其他为符号为已知数
(r1-s)^2 / (r2-s)^2 - (h1-h2)^2 / b^2 - 1 = 0
(r3-s)^2 / (r2-s)^2 - (h3-h2)^2 / b^2 - 1 = 0
[b*(r3-s) / (r2-s)^2 ] / [(r3-s)^2 / (r2-s)^2-1)]^0.5 - k = 0

以往是用的 matlab 的 solve 函数来求解方程组的符号解，现在希望脱离 matlab 环境，程序能够独立运行。
最好能直接使用 C#能够调用的数学库，当前搜集到的方法有：

1.利用 Matlab 二次开发，需要电脑上安装有 Matlab 或者 mablab 环境，不方便
2.C#有第三方的 mathnet 数学库，但是只能用来解线性方程，没有非线性方程的功能
3.C++的 gsl 库，貌似是利用牛顿迭代法求的根，但只能得到一个解，并且需要事先指定初始值
4.python 的相关数学库，没有试过

请大家提供下思路，还有有什么库或者方法，最好是能像 matlab 那样计算出所有的符号解或者多个数值解来，再次感谢。

[shiji]

我记得 MATLAB 有一个模组可以把 MATLAB 代码转换成 C 语言代码，你可以深入研究一下

[kindlebeta]

@shiji 是可以脱离 matlab 环境而独立运行吗？我找找看， 3q

[hexasnake]

@shiji 曾经我也这么尝试过，结果失败了。大概是 4 年前的事情。

[shiji]

@kindlebeta http://www.mathworks.com/products/matlab-coder/ http://www.mathworks.com/products/matlab-coder/features.html#generating-code-from-matlab


@hexasnake 可能那个时候 Matlab 这个功能还不成熟吧

[kindlebeta]

@shiji 难道要收费，不过好像有试用，试试看，再次感谢您

[ericls]

numpy sympy statsmodel

[ericls]

符号运算 sympy

[theoractice]

@shiji 你乐观了。问题在于并非所有函数都支持 codegen 。。。虽然 matlab 每个新版本都会把一些老函数开源，但我猜 mathworks 是永远不会让 solve 这种核心 API 支持 codegen 的（至少 R2015b 仍然不行），不然它还怎么赚钱。

另外既然是工程软件， LZ 干脆用 matlab 多算几种可能的几种符号解，预置进去给用户选好了。 sympy 估计是现有的最好方案，但说实话效果很差。变变花样它就未必解得出来。

[shiji]

@theoractice 我也没用过，只记得有这么个东西。


@kindlebeta 能支持转换 C 的函数列表在此， http://www.mathworks.com/help/coder/ug/functions-supported-for-code-generation--alphabetical-list.html?refresh=true

[theoractice]

其实我还真试了一下 sympy ，半小时了还没出结果哈哈哈哈
import sympy as sp

s, b, h2 = symbols('s b h2')
r1, r2 ,r3 = symbols('r1 r2 r3')
k, h1, h3 = symbols('k h1 h3')

eq1 = (r1-s)**2 / (r2-s)**2 - (h1-h2)**2 / b**2 - 1
eq2 = (r3-s)**2 / (r2-s)**2 - (h3-h2)**2 / b**2 - 1
eq3 = (b*(r3-s) / (r2-s)**2) / ((r3-s)**2 / (r2-s)**2-1)**0.5 - k

sp.solvers.solve((eq1, eq2, eq3), (s, b, h2))

[northisland]

问题是 solve [multivariant] non-linear equations

搬运 http://stackoverflow.com/questions/4233964/what-good-libraries-are-there-for-solving-a-system-of-non-linear-equations-in-c 中所说：

sundials packages
https://computation.llnl.gov/casc/sundials/description/description.html#descr_kinsol
https://computation.llnl.gov/casc/sundials/documentation/kin_examples.pdf

NLEQ
http://elib.zib.de/pub/elib/codelib/NewtonLib/index.html

继续围观

[bigtan]

MATLAB 的代码编译成 C/C#/Java 的库是可行的，我实践过 C#和 java 。

[theoractice]

结果出来了，我就不贴了，贴出来占两个半屏幕。 LZ 你确定真的要这么干么。。。

[jakiepaper]

@kindlebeta 有 r1, r2 ,r3 ， k, h1, h3 的具体数值吗？我可以用 sundials 帮你试试。我觉得非线性的方程组都很难解的，特别是在没有一个好的 initial guess 的情况下。

[kindlebeta]

@theoractice 真的能算出符号解吗，您能否把结果或者计算的方法发一下吗？我用 mathmetica 怎么都算不出来。佩服佩服

[kindlebeta]

r1=20.25 r2=19.5 r3=25.54 h1=88.8 k=-3.95
您看看这组数据可以算不
下面是解
b s h2
实数 虚数 实数 虚数 实数 虚数
8.52413E-14 115.7121688 27.82349676 0 136.8022229 0
-167.6813172 0 -137.9741991 0 72.41516435 0
-7.58524E-14 -42.62336039 23.8914588 0 64.97646733 0
-9.71682E-14 -133.8518931 24.74488104 0 157.7753132 0

[kindlebeta]

@jakiepaper 刚才少了一个 h3 h3=25.5

[theoractice]

matlab 算的啊，你难道没有试么？还是我算错了。

>> [s b h2]=solve('(r1-s)^2 / (r2-s)^2 - (h1-h2)^2 / b^2 - 1','(r3-s)^2 / (r2-s)^2 - (h3-h2)^2 / b^2 - 1','(b*(r3-s) / (r2-s)^2) / ((r3-s)^2 / (r2-s)^2-1)^0.5 - k','s','b','h2')

[theoractice]

你那个是数值解吧。用你的数据验证了下，结果有偏差但不大。你用 matlab 跑一下我上面那行，等个十几分钟就能出结果。

[cbsw]

http://www.sagemath.org/

[kindlebeta]

@theoractice 当时算了好长时间都没反应，还以为不能算呢，哈哈

[theoractice]

@kindlebeta 这种科学计算一俩小时算不出来太正常了。只要没出错，等一天也得等。

[kindlebeta]

@theoractice 我用你上面的算了一下，却算不出来，提示的却是
Warning: Explicit solution could not be found.
> In solve at 81
s =
[ empty sym ]
b =
[]
h2 =
[]
>>
难道是因为版本的问题吗？我的是 R2010a

[theoractice]

那必须是版本问题。我来刷个屏
>> [s b h2]=solve('(r1-s)^2 / (r2-s)^2 - (h1-h2)^2 / b^2 - 1','(r3-s)^2 / (r2-s)^2 - (h3-h2)^2 / b^2 - 1','(b*(r3-s) / (r2-s)^2) / ((r3-s)^2 / (r2-s)^2-1)^0.5 - k','s','b','h2')

s =

(h1^2*r3 - 2.0*h1*h3*r3 - 2.0*h1*k*r2^2 + 2.0*h1*k*r3^2 + h3^2*r3 - 1.0*h3*k*r1^2 + 2.0*h3*k*r1*r3 + 2.0*h3*k*r2^2 - 3.0*h3*k*r3^2 + 2.0*k^2*r1*r2^2 - 2.0*k^2*r1*r3^2 - 2.0*k^2*r2^2*r3 + 2.0*k^2*r3^3)/(h1^2 - 2.0*h1*h3 + h3^2 + 4.0*k^2*r3^2 + 4.0*k^2*r1*r2 - 4.0*k^2*r1*r3 - 4.0*k^2*r2*r3 - 4.0*h1*k*r2 + 4.0*h1*k*r3 + 4.0*h3*k*r2 - 4.0*h3*k*r3) + ((k*r1^2 - 2.0*k*r1*r3 + k*r3^2)*(h1*r2^2 + h1*r3^2 + h3*r1^2 - 1.0*h3*r2^2 + r2*(-1.0*(r1 - 1.0*r2)*(h1 - 1.0*h3 - 1.0*k*r1 + k*r3)*(h1*r1 + h1*r2 - 2.0*h1*r3 - 1.0*h3*r1 - 1.0*h3*r2 + 2.0*h3*r3 + k*r1^2 - 1.0*k*r1*r2 - 1.0*k*r1*r3 + k*r2*r3))^(1/2) - 1.0*r3*(-1.0*(r1 - 1.0*r2)*(h1 - 1.0*h3 - 1.0*k*r1 + k*r3)*(h1*r1 + h1*r2 - 2.0*h1*r3 - 1.0*h3*r1 - 1.0*h3*r2 + 2.0*h3*r3 + k*r1^2 - 1.0*k*r1*r2 - 1.0*k*r1*r3 + k*r2*r3))^(1/2) - 1.0*k*r1*r2^2 + k*r1^2*r2 + k*r1*r3^2 - 1.0*k*r1^2*r3 - 1.0*k*r2*r3^2 + k*r2^2*r3 - 2.0*h1*r2*r3 - 2.0*h3*r1*r3 + 2.0*h3*r2*r3))/((r1^2 - 2.0*r1*r3 + r3^2)*(h1^2 - 2.0*h1*h3 + h3^2 + 4.0*k^2*r3^2 + 4.0*k^2*r1*r2 - 4.0*k^2*r1*r3 - 4.0*k^2*r2*r3 - 4.0*h1*k*r2 + 4.0*h1*k*r3 + 4.0*h3*k*r2 - 4.0*h3*k*r3))
(h1^2*r3 - 2.0*h1*h3*r3 - 2.0*h1*k*r2^2 + 2.0*h1*k*r3^2 + h3^2*r3 - 1.0*h3*k*r1^2 + 2.0*h3*k*r1*r3 + 2.0*h3*k*r2^2 - 3.0*h3*k*r3^2 + 2.0*k^2*r1*r2^2 - 2.0*k^2*r1*r3^2 - 2.0*k^2*r2^2*r3 + 2.0*k^2*r3^3)/(h1^2 - 2.0*h1*h3 + h3^2 + 4.0*k^2*r3^2 + 4.0*k^2*r1*r2 - 4.0*k^2*r1*r3 - 4.0*k^2*r2*r3 - 4.0*h1*k*r2 + 4.0*h1*k*r3 + 4.0*h3*k*r2 - 4.0*h3*k*r3) + ((k*r1^2 - 2.0*k*r1*r3 + k*r3^2)*(h1*r2^2 + h1*r3^2 + h3*r1^2 - 1.0*h3*r2^2 - 1.0*r2*(-1.0*(r1 - 1.0*r2)*(h1 - 1.0*h3 - 1.0*k*r1 + k*r3)*(h1*r1 + h1*r2 - 2.0*h1*r3 - 1.0*h3*r1 - 1.0*h3*r2 + 2.0*h3*r3 + k*r1^2 - 1.0*k*r1*r2 - 1.0*k*r1*r3 + k*r2*r3))^(1/2) + r3*(-1.0*(r1 - 1.0*r2)*(h1 - 1.0*h3 - 1.0*k*r1 + k*r3)*(h1*r1 + h1*r2 - 2.0*h1*r3 - 1.0*h3*r1 - 1.0*h3*r2 + 2.0*h3*r3 + k*r1^2 - 1.0*k*r1*r2 - 1.0*k*r1*r3 + k*r2*r3))^(1/2) - 1.0*k*r1*r2^2 + k*r1^2*r2 + k*r1*r3^2 - 1.0*k*r1^2*r3 - 1.0*k*r2*r3^2 + k*r2^2*r3 - 2.0*h1*r2*r3 - 2.0*h3*r1*r3 + 2.0*h3*r2*r3))/((r1^2 - 2.0*r1*r3 + r3^2)*(h1^2 - 2.0*h1*h3 + h3^2 + 4.0*k^2*r3^2 + 4.0*k^2*r1*r2 - 4.0*k^2*r1*r3 - 4.0*k^2*r2*r3 - 4.0*h1*k*r2 + 4.0*h1*k*r3 + 4.0*h3*k*r2 - 4.0*h3*k*r3))
(h1^2*r3 - 2.0*h1*h3*r3 + 2.0*h1*k*r2^2 - 2.0*h1*k*r3^2 + h3^2*r3 + h3*k*r1^2 - 2.0*h3*k*r1*r3 - 2.0*h3*k*r2^2 + 3.0*h3*k*r3^2 + 2.0*k^2*r1*r2^2 - 2.0*k^2*r1*r3^2 - 2.0*k^2*r2^2*r3 + 2.0*k^2*r3^3)/(h1^2 - 2.0*h1*h3 + h3^2 + 4.0*k^2*r3^2 + 4.0*k^2*r1*r2 - 4.0*k^2*r1*r3 - 4.0*k^2*r2*r3 + 4.0*h1*k*r2 - 4.0*h1*k*r3 - 4.0*h3*k*r2 + 4.0*h3*k*r3) - (1.0*(k*r1^2 - 2.0*k*r1*r3 + k*r3^2)*(h1*r2^2 + h1*r3^2 + h3*r1^2 - 1.0*h3*r2^2 + r2*(-1.0*(r1 - 1.0*r2)*(h1 - 1.0*h3 + k*r1 - 1.0*k*r3)*(h1*r1 + h1*r2 - 2.0*h1*r3 - 1.0*h3*r1 - 1.0*h3*r2 + 2.0*h3*r3 - 1.0*k*r1^2 + k*r1*r2 + k*r1*r3 - 1.0*k*r2*r3))^(1/2) - 1.0*r3*(-1.0*(r1 - 1.0*r2)*(h1 - 1.0*h3 + k*r1 - 1.0*k*r3)*(h1*r1 + h1*r2 - 2.0*h1*r3 - 1.0*h3*r1 - 1.0*h3*r2 + 2.0*h3*r3 - 1.0*k*r1^2 + k*r1*r2 + k*r1*r3 - 1.0*k*r2*r3))^(1/2) + k*r1*r2^2 - 1.0*k*r1^2*r2 - 1.0*k*r1*r3^2 + k*r1^2*r3 + k*r2*r3^2 - 1.0*k*r2^2*r3 - 2.0*h1*r2*r3 - 2.0*h3*r1*r3 + 2.0*h3*r2*r3))/((r1^2 - 2.0*r1*r3 + r3^2)*(h1^2 - 2.0*h1*h3 + h3^2 + 4.0*k^2*r3^2 + 4.0*k^2*r1*r2 - 4.0*k^2*r1*r3 - 4.0*k^2*r2*r3 + 4.0*h1*k*r2 - 4.0*h1*k*r3 - 4.0*h3*k*r2 + 4.0*h3*k*r3))
(h1^2*r3 - 2.0*h1*h3*r3 + 2.0*h1*k*r2^2 - 2.0*h1*k*r3^2 + h3^2*r3 + h3*k*r1^2 - 2.0*h3*k*r1*r3 - 2.0*h3*k*r2^2 + 3.0*h3*k*r3^2 + 2.0*k^2*r1*r2^2 - 2.0*k^2*r1*r3^2 - 2.0*k^2*r2^2*r3 + 2.0*k^2*r3^3)/(h1^2 - 2.0*h1*h3 + h3^2 + 4.0*k^2*r3^2 + 4.0*k^2*r1*r2 - 4.0*k^2*r1*r3 - 4.0*k^2*r2*r3 + 4.0*h1*k*r2 - 4.0*h1*k*r3 - 4.0*h3*k*r2 + 4.0*h3*k*r3) - (1.0*(k*r1^2 - 2.0*k*r1*r3 + k*r3^2)*(h1*r2^2 + h1*r3^2 + h3*r1^2 - 1.0*h3*r2^2 - 1.0*r2*(-1.0*(r1 - 1.0*r2)*(h1 - 1.0*h3 + k*r1 - 1.0*k*r3)*(h1*r1 + h1*r2 - 2.0*h1*r3 - 1.0*h3*r1 - 1.0*h3*r2 + 2.0*h3*r3 - 1.0*k*r1^2 + k*r1*r2 + k*r1*r3 - 1.0*k*r2*r3))^(1/2) + r3*(-1.0*(r1 - 1.0*r2)*(h1 - 1.0*h3 + k*r1 - 1.0*k*r3)*(h1*r1 + h1*r2 - 2.0*h1*r3 - 1.0*h3*r1 - 1.0*h3*r2 + 2.0*h3*r3 - 1.0*k*r1^2 + k*r1*r2 + k*r1*r3 - 1.0*k*r2*r3))^(1/2) + k*r1*r2^2 - 1.0*k*r1^2*r2 - 1.0*k*r1*r3^2 + k*r1^2*r3 + k*r2*r3^2 - 1.0*k*r2^2*r3 - 2.0*h1*r2*r3 - 2.0*h3*r1*r3 + 2.0*h3*r2*r3))/((r1^2 - 2.0*r1*r3 + r3^2)*(h1^2 - 2.0*h1*h3 + h3^2 + 4.0*k^2*r3^2 + 4.0*k^2*r1*r2 - 4.0*k^2*r1*r3 - 4.0*k^2*r2*r3 + 4.0*h1*k*r2 - 4.0*h1*k*r3 - 4.0*h3*k*r2 + 4.0*h3*k*r3))

[theoractice]

哎， b 实在太长这里贴不下。还是不污染环境了。
http://pastebin.com/ddDcu0N7
LZ 去看吧。

[facat]

用牛顿法吧，自己写很容易的

[jakiepaper]

@kindlebeta 我不会。。。（哭）没解过复数的问题，不会。