gpt4 book ai didi

math - 在 Mathematica 中求解非线性方程组

转载 作者:行者123 更新时间:2023-12-01 23:33:59 26 4
gpt4 key购买 nike

我正在尝试解决 Mathemtica 中的非线性方程组。
我尝试了 Solve 和 NSolve,我也尝试定义 a_{ij} 和 b_{ij} 和 m33=1 数值来简化方程,但是 Mathematica 似乎工作太久或者我做错了什么。在 Mathematica 我只是想找到解决方案,但我还需要一些 c/c++ 库来在我的代码中执行此操作。

“运算符”中的主要方程:

M[A[(x,y)]]=B[M[(x,y)]]

其中“运算符”是透视变换:
u= (m13 + m11*x + m12*y)/(m33 + m31*x + m32*y); 

v= (m23 + m21*x +m22*y)/(m33 + m31*x + m32*y);

我在 Mathematica 中的输入:
Solve[(b13 + (b11 (m13 + m11 x1 + m12 y1))/(m33 + m31 x1 + 
m32 y1) + (b12 (m23 + m21 x1 + m22 y1))/(m33 + m31 x1 +
m32 y1))/(b33 + (b31 (m13 + m11 x1 + m12 y1))/(m33 + m31 x1 +
m32 y1) + (b32 (m23 + m21 x1 + m22 y1))/(m33 + m31 x1 +
m32 y1)) == (m13 + (m11 (a13 + a11 x1 + a12 y1))/(a33 +
a31 x1 + a32 y1) + (m12 (a23 + a21 x1 + a22 y1))/(a33 +
a31 x1 + a32 y1))/(m33 + (m31 (a13 + a11 x1 + a12 y1))/(a33 +
a31 x1 + a32 y1) + (m32 (a23 + a21 x1 + a22 y1))/(a33 +
a31 x1 +
a32 y1)) && (b23 + (b21 (m13 + m11 x1 + m12 y1))/(m33 +
m31 x1 + m32 y1) + (b22 (m23 + m21 x1 + m22 y1))/(m33 +
m31 x1 + m32 y1))/(b33 + (b31 (m13 + m11 x1 + m12 y1))/(m33 +
m31 x1 + m32 y1) + (b32 (m23 + m21 x1 + m22 y1))/(m33 +
m31 x1 +
m32 y1)) == (m23 + (m21 (a13 + a11 x1 + a12 y1))/(a33 +
a31 x1 + a32 y1) + (m22 (a23 + a21 x1 + a22 y1))/(a33 +
a31 x1 + a32 y1))/(m33 + (m31 (a13 + a11 x1 + a12 y1))/(a33 +
a31 x1 + a32 y1) + (m32 (a23 + a21 x1 + a22 y1))/(a33 +
a31 x1 +
a32 y1)) && (b13 + (b11 (m13 + m11 x2 + m12 y2))/(m33 +
m31 x2 + m32 y2) + (b12 (m23 + m21 x2 + m22 y2))/(m33 +
m31 x2 + m32 y2))/(b33 + (b31 (m13 + m11 x2 + m12 y2))/(m33 +
m31 x2 + m32 y2) + (b32 (m23 + m21 x2 + m22 y2))/(m33 +
m31 x2 +
m32 y2)) == (m13 + (m11 (a13 + a11 x2 + a12 y2))/(a33 +
a31 x2 + a32 y2) + (m12 (a23 + a21 x2 + a22 y2))/(a33 +
a31 x2 + a32 y2))/(m33 + (m31 (a13 + a11 x2 + a12 y2))/(a33 +
a31 x2 + a32 y2) + (m32 (a23 + a21 x2 + a22 y2))/(a33 +
a31 x2 +
a32 y2)) && (b23 + (b21 (m13 + m11 x2 + m12 y2))/(m33 +
m31 x2 + m32 y2) + (b22 (m23 + m21 x2 + m22 y2))/(m33 +
m31 x2 + m32 y2))/(b33 + (b31 (m13 + m11 x2 + m12 y2))/(m33 +
m31 x2 + m32 y2) + (b32 (m23 + m21 x2 + m22 y2))/(m33 +
m31 x2 +
m32 y2)) == (m23 + (m21 (a13 + a11 x2 + a12 y2))/(a33 +
a31 x2 + a32 y2) + (m22 (a23 + a21 x2 + a22 y2))/(a33 +
a31 x2 + a32 y2))/(m33 + (m31 (a13 + a11 x2 + a12 y2))/(a33 +
a31 x2 + a32 y2) + (m32 (a23 + a21 x2 + a22 y2))/(a33 +
a31 x2 +
a32 y2)) && (b13 + (b11 (m13 + m11 x3 + m12 y3))/(m33 +
m31 x3 + m32 y3) + (b12 (m23 + m21 x3 + m22 y3))/(m33 +
m31 x3 + m32 y3))/(b33 + (b31 (m13 + m11 x3 + m12 y3))/(m33 +
m31 x3 + m32 y3) + (b32 (m23 + m21 x3 + m22 y3))/(m33 +
m31 x3 +
m32 y3)) == (m13 + (m11 (a13 + a11 x3 + a12 y3))/(a33 +
a31 x3 + a32 y3) + (m12 (a23 + a21 x3 + a22 y3))/(a33 +
a31 x3 + a32 y3))/(m33 + (m31 (a13 + a11 x3 + a12 y3))/(a33 +
a31 x3 + a32 y3) + (m32 (a23 + a21 x3 + a22 y3))/(a33 +
a31 x3 +
a32 y3)) && (b23 + (b21 (m13 + m11 x3 + m12 y3))/(m33 +
m31 x3 + m32 y3) + (b22 (m23 + m21 x3 + m22 y3))/(m33 +
m31 x3 + m32 y3))/(b33 + (b31 (m13 + m11 x3 + m12 y3))/(m33 +
m31 x3 + m32 y3) + (b32 (m23 + m21 x3 + m22 y3))/(m33 +
m31 x3 +
m32 y3)) == (m23 + (m21 (a13 + a11 x3 + a12 y3))/(a33 +
a31 x3 + a32 y3) + (m22 (a23 + a21 x3 + a22 y3))/(a33 +
a31 x3 + a32 y3))/(m33 + (m31 (a13 + a11 x3 + a12 y3))/(a33 +
a31 x3 + a32 y3) + (m32 (a23 + a21 x3 + a22 y3))/(a33 +
a31 x3 +
a32 y3)) && (b13 + (b11 (m13 + m11 x4 + m12 y4))/(m33 +
m31 x4 + m32 y4) + (b12 (m23 + m21 x4 + m22 y4))/(m33 +
m31 x4 + m32 y4))/(b33 + (b31 (m13 + m11 x4 + m12 y4))/(m33 +
m31 x4 + m32 y4) + (b32 (m23 + m21 x4 + m22 y4))/(m33 +
m31 x4 +
m32 y4)) == (m13 + (m11 (a13 + a11 x4 + a12 y4))/(a33 +
a31 x4 + a32 y4) + (m12 (a23 + a21 x4 + a22 y4))/(a33 +
a31 x4 + a32 y4))/(m33 + (m31 (a13 + a11 x4 + a12 y4))/(a33 +
a31 x4 + a32 y4) + (m32 (a23 + a21 x4 + a22 y4))/(a33 +
a31 x4 +
a32 y4)) && (b23 + (b21 (m13 + m11 x4 + m12 y4))/(m33 +
m31 x4 + m32 y4) + (b22 (m23 + m21 x4 + m22 y4))/(m33 +
m31 x4 + m32 y4))/(b33 + (b31 (m13 + m11 x4 + m12 y4))/(m33 +
m31 x4 + m32 y4) + (b32 (m23 + m21 x4 + m22 y4))/(m33 +
m31 x4 +
m32 y4)) == (m23 + (m21 (a13 + a11 x4 + a12 y4))/(a33 +
a31 x4 + a32 y4) + (m22 (a23 + a21 x4 + a22 y4))/(a33 +
a31 x4 + a32 y4))/(m33 + (m31 (a13 + a11 x4 + a12 y4))/(a33 +
a31 x4 + a32 y4) + (m32 (a23 + a21 x4 + a22 y4))/(a33 +
a31 x4 + a32 y4)) && m33 == 1, {m11, m12, m13, m21, m22, m23,
m31, m32}]

最佳答案

为了有任何机会,您将需要实际数字代替参数变量。否则系统太大而无法处理。

为了说明这一点,我创建了多项式(忽略分母消失的场景),然后对参数进行了随机数字替换。我本可以通过替换消除 m33,但选择在系统中保留 m33-1==0。这样就不需要对任何一个方程进行特殊处理。为了提高效率,人们可能会考虑对变量线性的方程子集进行这种消除。

In[40]:= exprs = Apply[Subtract, eqns, {1}];
e2 = Together[exprs];
polys = Numerator[e2];

In[62]:= allvars = Variables[polys];
vars = {m11, m12, m13, m21, m22, m23, m31, m32, m33};
params = Complement[allvars, vars]

Out[64]= {a11, a12, a13, a21, a22, a23, a31, a32, a33, b11, b12, b13, \
b21, b22, b23, b31, b32, b33, x1, x2, x3, x4, y1, y2, y3, y4}

In[69]:= SeedRandom[11111];
substitutions =
Thread[params -> RandomInteger[{-1000, 1000}, Length[params]]];

In[71]:= numpolys = polys /. substitutions;

NSolve 认为该系统实际上是欠定的,也就是说,您的方程具有冗余(代数相关,更专业地说)。所以它与一个伪随机超平面相交,然后得到一个有限解集。
In[73]:= Timing[solns = NSolve[numpolys == 0, vars];]

During evaluation of In[73]:= NSolve::infsolns: Infinite solution set has dimension at least 1. Returning intersection of solutions with (107814 m11)/118505-(177066 m12)/118505-(164294 m13)/118505+(32943 m21)/23701+(186238 m22)/118505-(126102 m23)/118505-(178233 m31)/118505-(185338 m32)/118505+(141088 m33)/118505 == 1.

Out[73]= {357.420000, Null}

以下是这种情况下的解决方案。
In[74]:= solns // N

Out[74]= {{m11 -> -2.22241, m12 -> 0., m13 -> -2.41203,
m21 -> -0.539924, m22 -> 2.33146*10^-172, m23 -> -0.585993,
m31 -> 0.921382, m32 -> -4.05984*10^-172,
m33 -> 1.}, {m11 -> -2.22241, m12 -> 0., m13 -> -2.41203,
m21 -> -0.539924, m22 -> 2.33146*10^-172, m23 -> -0.585993,
m31 -> 0.921382, m32 -> -4.05984*10^-172,
m33 -> 1.}, {m11 -> -2.22241, m12 -> 0., m13 -> -2.41203,
m21 -> -0.539924, m22 -> 2.33146*10^-172, m23 -> -0.585993,
m31 -> 0.921382, m32 -> -4.05984*10^-172,
m33 -> 1.}, {m11 -> -0.029351, m12 -> 0., m13 -> 0.409304,
m21 -> 0.0182228, m22 -> -2.19075*10^-169, m23 -> -0.25412,
m31 -> -0.0717095, m32 -> 2.05529*10^-169,
m33 -> 1.}, {m11 -> -0.029351, m12 -> 0., m13 -> 0.409304,
m21 -> 0.0182228, m22 -> -2.19075*10^-169, m23 -> -0.25412,
m31 -> -0.0717095, m32 -> 2.05529*10^-169,
m33 -> 1.}, {m11 -> -0.029351, m12 -> 0., m13 -> 0.409304,
m21 -> 0.0182228, m22 -> -2.19075*10^-169, m23 -> -0.25412,
m31 -> -0.0717095, m32 -> 2.05529*10^-169,
m33 -> 1.}, {m11 -> 0.541883, m12 -> 0., m13 -> -0.123031,
m21 -> -4.58369, m22 -> -5.60174*10^-170, m23 -> 1.0407,
m31 -> -4.40445, m32 -> -5.32622*10^-170,
m33 -> 1.}, {m11 -> 0.541883, m12 -> 0., m13 -> -0.123031,
m21 -> -4.58369, m22 -> -5.60174*10^-170, m23 -> 1.0407,
m31 -> -4.40445, m32 -> -5.32622*10^-170,
m33 -> 1.}, {m11 -> 0.541883, m12 -> 0., m13 -> -0.123031,
m21 -> -4.58369, m22 -> -5.60174*10^-170, m23 -> 1.0407,
m31 -> -4.40445, m32 -> -5.32622*10^-170, m33 -> 1.}}

我们检查它们是否满足具有非常小的数字残差的原始表达式。
In[76]:= exprs /. substitutions /. solns

Out[76]= {{-4.44089*10^-16, 1.11022*10^-16, -8.88178*10^-16, 0.,
0., -1.11022*10^-16, -4.44089*10^-16, 1.11022*10^-16,
0.}, {-4.44089*10^-16, 1.11022*10^-16, -8.88178*10^-16, 0.,
0., -1.11022*10^-16, -4.44089*10^-16, 1.11022*10^-16,
0.}, {-4.44089*10^-16, 1.11022*10^-16, -8.88178*10^-16, 0.,
0., -1.11022*10^-16, -4.44089*10^-16, 1.11022*10^-16,
0.}, {-2.22045*10^-16, 1.11022*10^-16, -1.66533*10^-16,
1.11022*10^-16, -5.55112*10^-17, 1.11022*10^-16, -1.11022*10^-16,
5.55112*10^-17, 0.}, {-2.22045*10^-16,
1.11022*10^-16, -1.66533*10^-16, 1.11022*10^-16, -5.55112*10^-17,
1.11022*10^-16, -1.11022*10^-16, 5.55112*10^-17,
0.}, {-2.22045*10^-16, 1.11022*10^-16, -1.66533*10^-16,
1.11022*10^-16, -5.55112*10^-17, 1.11022*10^-16, -1.11022*10^-16,
5.55112*10^-17, 0.}, {2.34535*10^-15, -1.11022*10^-15,
1.11022*10^-16, 2.22045*10^-16, 1.31839*10^-15, -1.11022*10^-15,
1.249*10^-15, -6.66134*10^-16,
0.}, {2.34535*10^-15, -1.11022*10^-15, 1.11022*10^-16,
2.22045*10^-16, 1.31839*10^-15, -1.11022*10^-15,
1.249*10^-15, -6.66134*10^-16,
0.}, {2.34535*10^-15, -1.11022*10^-15, 1.11022*10^-16,
2.22045*10^-16, 1.31839*10^-15, -1.11022*10^-15,
1.249*10^-15, -6.66134*10^-16, 0.}}

关于math - 在 Mathematica 中求解非线性方程组,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/11735078/

26 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com