注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Mr.Right

不顾一切的去想,于是我们有了梦想。脚踏实地的去做,于是梦想成了现实。

 
 
 

日志

 
 
关于我

人生一年又一年,只要每年都有所积累,有所成长,都有那么一次自己认为满意的花开时刻就好。即使一时不顺,也要敞开胸怀。生命的荣枯并不是简单的重复,一时的得失不是成败的尺度。花开不是荣耀,而是一个美丽的结束,花谢也不是耻辱,而是一个低调的开始。

网易考拉推荐

mathematica中Resultant使用范例  

2012-04-04 15:16:23|  分类: 学习 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

缘起:我以前不知道求两个多项式关于某个变量的Resultant可以消去该变量

Tips & Tricks:  ① %n 或 Out[n] 返回第 n 个输出的值 ② Ctrl+K是参数补齐,类似matlab中的Tab键③Ctrl+Shift+B匹配括号,光标放在括号左边

1)分析为什么下面语句执行后aZeros返回如下结果

poly1[z_, a_] = Product[z - i a, {i, 1, 7, 2}] + 2

poly2[z_, a_] = Product[z - i a, {i, 0, 6, 2}] + 3

Resultant[poly1[z, a], poly2[z, a], z]   (*消去z后关于a的式子*)

Resultant[poly1[z, a], poly2[z, a], a]

aZeros = Chop[Last /@ (List @@ NRoots[%% == 0, a])]   (* 求解关于a的表达式的数值解 *)

{-0.697333, -0.682883, -0.558864 - 0.558864 I, -0.558864 +
  0.558864 I, -0.0874929 - 0.0874929 I, -0.0874929 + 0.0874929 I,
 0. - 0.682883 I, 0. + 0.682883 I, 0. - 0.697333 I, 0. + 0.697333 I,
 0.0874929 - 0.0874929 I, 0.0874929 + 0.0874929 I,
 0.558864 - 0.558864 I, 0.558864 + 0.558864 I, 0.682883, 0.697333}

  1. 执行 NRoots[%3 == 0, a]
    返回 a == -0.697333 || a == -0.682883 || a == -0.558864 - 0.558864 I ||
     a == -0.558864 + 0.558864 I || a == -0.0874929 - 0.0874929 I ||
     a == -0.0874929 + 0.0874929 I || a == 0. - 0.682883 I ||
     a == 0. + 0.682883 I || a == 0. - 0.697333 I ||
     a == 0. + 0.697333 I || a == 0.0874929 - 0.0874929 I ||
     a == 0.0874929 + 0.0874929 I || a == 0.558864 - 0.558864 I ||
     a == 0.558864 + 0.558864 I || a == 0.682883 || a == 0.697333
  2. 执行List @@ NRoots[%% == 0, a]  (*原来是去掉中间的 || 算符,并合并成为集合*)
    返回 {a == -0.697333, a == -0.682883, a == -0.558864 - 0.558864 I,
     a == -0.558864 + 0.558864 I, a == -0.0874929 - 0.0874929 I,
     a == -0.0874929 + 0.0874929 I, a == 0. - 0.682883 I,
     a == 0. + 0.682883 I, a == 0. - 0.697333 I, a == 0. + 0.697333 I,
     a == 0.0874929 - 0.0874929 I, a == 0.0874929 + 0.0874929 I,
     a == 0.558864 - 0.558864 I, a == 0.558864 + 0.558864 I,
     a == 0.682883, a == 0.697333}
  3. 执行 Last /@ (List @@ NRoots[%3 == 0, a])  (* 看到Last的妙用了吗?取RHS*)
    返回 {-0.697333, -0.682883, -0.558864 - 0.558864 I, -0.558864 +
      0.558864 I, -0.0874929 - 0.0874929 I, -0.0874929 + 0.0874929 I,
     0. - 0.682883 I, 0. + 0.682883 I, 0. - 0.697333 I, 0. + 0.697333 I,
     0.0874929 - 0.0874929 I, 0.0874929 + 0.0874929 I,
     0.558864 - 0.558864 I, 0.558864 + 0.558864 I, 0.682883, 0.697333}
  4. 执行 Chop[Last /@ (List @@ NRoots[%% == 0, a])]
    返回的结果与不加Chop是一样的,Chop[ ]函数的作用是将近似为0的函数变为0

2) 结论是aZeros = Chop[Last /@ (List @@ NRoots[%% == 0, a])]的核心部分就是NRoots[%% == 0, a],其他部分都是修饰该结果的

3) 求关于变量z的Resultant的作用就是联立求解以下两个等式消去z,解出关于a的表达式等于零的根。下面是如何联立两个方程

执行

Sort @ Table[{a -> aZeros[[i]] , z -> Flatten[ (* 解集合中是{a,z}对, a取自于aZeros数组z原来就是后面这一串的计算结果*)  (* Sort @ 是 Sort[...] 的等价写法*)
Outer[
If[(* which are the same? *)
Chop[#1 - #2] === 0, Chop[#1], {}] &,   (* 这里为Pure Function,#1为第一个参数 *)
(* the zeros of poly1 with respect to z *)
Last /@ List @@ NRoots[poly1[z, aZeros[[i]]] == 0, z],
(* the zeros of poly2 with respect to z *)
Last /@ List @@ NRoots[poly2[z, aZeros[[i]]] == 0, z]]] [[1]] },
{i, Length[aZeros]}]  (* 以前没看出来这里是Table的i的范围 *)

返回

{{a -> -0.697333, z -> -0.848507}, {a -> -0.682883,
  z -> -3.78763}, {a -> -0.558864 - 0.558864 I,
  z -> -1.88035 - 1.88035 I}, {a -> -0.558864 + 0.558864 I,
  z -> -1.88035 + 1.88035 I}, {a -> -0.0874929 - 0.0874929 I,
  z -> -1.21371 - 1.21371 I}, {a -> -0.0874929 + 0.0874929 I,
  z -> -1.21371 + 1.21371 I}, {a -> 0. - 0.682883 I,
  z -> 0. - 3.78763 I}, {a -> 0. + 0.682883 I,
  z -> 0. + 3.78763 I}, {a -> 0. - 0.697333 I,
  z -> 0. - 0.848507 I}, {a -> 0. + 0.697333 I,
  z -> 0. + 0.848507 I}, {a -> 0.0874929 - 0.0874929 I,
  z -> 1.21371 - 1.21371 I}, {a -> 0.0874929 + 0.0874929 I,
  z -> 1.21371 + 1.21371 I}, {a -> 0.558864 - 0.558864 I,
  z -> 1.88035 - 1.88035 I}, {a -> 0.558864 + 0.558864 I,
  z -> 1.88035 + 1.88035 I}, {a -> 0.682883,
  z -> 3.78763}, {a -> 0.697333, z -> 0.848507}}

  1. 执行 Table[NRoots[poly1[z, aZeros[[i]]] == 0, z], {i, Length[aZeros]}]
    返回 {z == -4.73016 || z == -3.8362 || z == -1.74247 || z == -0.848507,
     z == -4.61514 || z == -3.78763 || z == -1.67544 || z == -0.84793,
     z == -3.9667 - 3.9667 I || z == -2.59056 - 2.59056 I ||
      z == -1.88035 - 1.88035 I || z == -0.504217 - 0.504217 I,
     z == -3.9667 + 3.9667 I || z == -2.59056 + 2.59056 I ||
      z == -1.88035 + 1.88035 I || z == -0.504217 + 0.504217 I,
     z == -1.21371 - 1.21371 I || z == -1.1682 + 0.468255 I ||
      z == 0.468255 - 1.1682 I || z == 0.513767 + 0.513767 I,
     z == -1.21371 + 1.21371 I || z == -1.1682 - 0.468255 I ||
      z == 0.468255 + 1.1682 I || z == 0.513767 - 0.513767 I,
     z == 0. - 0.84793 I || z == 0. - 1.67544 I || z == 0. - 3.78763 I ||
      z == 0. - 4.61514 I,
     z == 0. + 0.84793 I || z == 0. + 1.67544 I || z == 0. + 3.78763 I ||
      z == 0. + 4.61514 I,
     z == 0. - 0.848507 I || z == 0. - 1.74247 I || z == 0. - 3.8362 I ||
      z == 0. - 4.73016 I,
     z == 0. + 0.848507 I || z == 0. + 1.74247 I || z == 0. + 3.8362 I ||
      z == 0. + 4.73016 I,
     z == -0.513767 + 0.513767 I || z == -0.468255 - 1.1682 I ||
      z == 1.1682 + 0.468255 I || z == 1.21371 - 1.21371 I,
     z == -0.513767 - 0.513767 I || z == -0.468255 + 1.1682 I ||
      z == 1.1682 - 0.468255 I || z == 1.21371 + 1.21371 I,
     z == 0.504217 - 0.504217 I || z == 1.88035 - 1.88035 I ||
      z == 2.59056 - 2.59056 I || z == 3.9667 - 3.9667 I,
     z == 0.504217 + 0.504217 I || z == 1.88035 + 1.88035 I ||
      z == 2.59056 + 2.59056 I || z == 3.9667 + 3.9667 I,
     z == 0.84793 || z == 1.67544 || z == 3.78763 || z == 4.61514,
     z == 0.848507 || z == 1.74247 || z == 3.8362 || z == 4.73016}
     
  2. 执行 Table[NRoots[poly2[z, aZeros[[i]]] == 0, z], {i, Length[aZeros]}]
    返回 {z == -3.91312 || z == -3.33549 || z == -0.848507 || z == -0.270884,
     z == -3.78763 || z == -3.32898 || z == -0.768318 || z == -0.309672,
     z == -3.43209 - 3.43209 I || z == -1.88035 - 1.88035 I ||
      z == -1.47284 - 1.47284 I || z == 0.0789058 + 0.0789058 I,
     z == -3.43209 + 3.43209 I || z == -1.88035 + 1.88035 I ||
      z == -1.47284 + 1.47284 I || z == 0.0789058 - 0.0789058 I,
     z == -1.21371 - 1.21371 I || z == -1.17258 + 0.647626 I ||
      z == 0.647626 - 1.17258 I || z == 0.688753 + 0.688753 I,
     z == -1.21371 + 1.21371 I || z == -1.17258 - 0.647626 I ||
      z == 0.647626 + 1.17258 I || z == 0.688753 - 0.688753 I,
     z == 0. - 0.309672 I || z == 0. - 0.768318 I || z == 0. - 3.32898 I ||
       z == 0. - 3.78763 I,
     z == 0. + 0.309672 I || z == 0. + 0.768318 I || z == 0. + 3.32898 I ||
       z == 0. + 3.78763 I,
     z == 0. - 0.270884 I || z == 0. - 0.848507 I || z == 0. - 3.33549 I ||
       z == 0. - 3.91312 I,
     z == 0. + 0.270884 I || z == 0. + 0.848507 I || z == 0. + 3.33549 I ||
       z == 0. + 3.91312 I,
     z == -0.688753 + 0.688753 I || z == -0.647626 - 1.17258 I ||
      z == 1.17258 + 0.647626 I || z == 1.21371 - 1.21371 I,
     z == -0.688753 - 0.688753 I || z == -0.647626 + 1.17258 I ||
      z == 1.17258 - 0.647626 I || z == 1.21371 + 1.21371 I,
     z == -0.0789058 + 0.0789058 I || z == 1.47284 - 1.47284 I ||
      z == 1.88035 - 1.88035 I || z == 3.43209 - 3.43209 I,
     z == -0.0789058 - 0.0789058 I || z == 1.47284 + 1.47284 I ||
      z == 1.88035 + 1.88035 I || z == 3.43209 + 3.43209 I,
     z == 0.309672 || z == 0.768318 || z == 3.32898 || z == 3.78763,
     z == 0.270884 || z == 0.848507 || z == 3.33549 || z == 3.91312}
  3. 执行List @@  Table[NRoots[poly1[z, aZeros[[i]]] == 0, z], {i, Length[aZeros]}]返回的结果与 Table[NRoots[poly1[z, aZeros[[i]]] == 0, z], {i, Length[aZeros]}] 完全相同,这说明给Table[NRoots[poly1[z, aZeros[[i]]] == 0, z], {i, Length[aZeros]}] Apply(即@@)一个List是重复的,因为Table[NRoots[poly1[z, aZeros[[i]]] == 0, z], {i, Length[aZeros]}]返回的结果本身就是一个集合。
  4. 执行 Last /@ List @@ Table[NRoots[poly1[z, aZeros[[i]]] == 0, z], {i, Length[aZeros]}]  或者执行 Last /@ Table[NRoots[poly1[z, aZeros[[i]]] == 0, z], {i, Length[aZeros]}],可以看出Last函数去掉了"||"
    返回 {z == -0.848507, z == -0.84793, z == -0.504217 - 0.504217 I,
     z == -0.504217 + 0.504217 I, z == 0.513767 + 0.513767 I,
     z == 0.513767 - 0.513767 I, z == 0. - 4.61514 I, z == 0. + 4.61514 I,
      z == 0. - 4.73016 I, z == 0. + 4.73016 I, z == 1.21371 - 1.21371 I,
     z == 1.21371 + 1.21371 I, z == 3.9667 - 3.9667 I,
     z == 3.9667 + 3.9667 I, z == 4.61514, z == 4.73016}
  5. 执行 g = If[(*which are the same?*)Chop[#1 - #2] === 0, Chop[#1], {}] &;
    g[2,3]返回{},g[2,2]返回 2 。
  6. 搞懂 Outer[f, list1, list2] 给出list1的广义外积,形成列表最底层元素的所有可能组合,并把它们作为 f 的自变量。比如执行 Outer[Times, {1, 2, 3, 4}, {a, b, c}] ,由于f是Times所以返回{{a, b, c}, {2 a, 2 b, 2 c}, {3 a, 3 b, 3 c}, {4 a, 4 b, 4 c}}
  7. 曾经困扰我的一点是没有看出来If[(*which are the same?*)Chop[#1 - #2] === 0, Chop[#1], {}] &原来是一个Pure Function, Outer[ If[(*which are the same?*)Chop[#1 - #2] === 0, Chop[#1], {}] &,
     Last /@ List @@ NRoots[poly1[z, aZeros[[i]]] == 0, z], 
    Last /@ List @@ NRoots[poly2[z, aZeros[[i]]] == 0, z] ]
  8. 执行 Table[Outer[If[Chop[#1 - #2] === 0, Chop[#1], {}] &,
      Last /@ List @@ NRoots[poly1[z, aZeros[[i]]] == 0, z],
      Last /@ List @@ NRoots[poly2[z, aZeros[[i]]] == 0, z]], {i,
      Length[aZeros]}]
    返回 {{{{}, {}, {}, {}}, {{}, {}, {}, {}}, {{}, {}, {}, {}}, {{}, {}, \
    -0.848507, {}}}, {{{}, {}, {}, {}}, {-3.78763, {}, {}, {}}, {{}, {}, \
    {}, {}}, {{}, {}, {}, {}}}, {{{}, {}, {}, {}}, {{}, {}, {}, {}}, {{}, \
    -1.88035 -
        1.88035 I, {}, {}}, {{}, {}, {}, {}}}, {{{}, {}, {}, {}}, {{}, \
    {}, {}, {}}, {{}, -1.88035 +
        1.88035 I, {}, {}}, {{}, {}, {}, {}}}, {{-1.21371 -
        1.21371 I, {}, {}, {}}, {{}, {}, {}, {}}, {{}, {}, {}, {}}, {{}, \
    {}, {}, {}}}, {{-1.21371 +
        1.21371 I, {}, {}, {}}, {{}, {}, {}, {}}, {{}, {}, {}, {}}, {{}, \
    {}, {}, {}}}, {{{}, {}, {}, {}}, {{}, {}, {}, {}}, {{}, {}, {},
       0. - 3.78763 I}, {{}, {}, {}, {}}}, {{{}, {}, {}, {}}, {{}, {}, \
    {}, {}}, {{}, {}, {}, 0. + 3.78763 I}, {{}, {}, {}, {}}}, {{{},
       0. - 0.848507 I, {}, {}}, {{}, {}, {}, {}}, {{}, {}, {}, {}}, {{}, \
    {}, {}, {}}}, {{{},
       0. + 0.848507 I, {}, {}}, {{}, {}, {}, {}}, {{}, {}, {}, {}}, {{}, \
    {}, {}, {}}}, {{{}, {}, {}, {}}, {{}, {}, {}, {}}, {{}, {}, {}, {}}, \
    {{}, {}, {},
       1.21371 -
        1.21371 I}}, {{{}, {}, {}, {}}, {{}, {}, {}, {}}, {{}, {}, {}, \
    {}}, {{}, {}, {}, 1.21371 + 1.21371 I}}, {{{}, {}, {}, {}}, {{}, {},
       1.88035 -
        1.88035 I, {}}, {{}, {}, {}, {}}, {{}, {}, {}, {}}}, {{{}, {}, \
    {}, {}}, {{}, {},
       1.88035 +
        1.88035 I, {}}, {{}, {}, {}, {}}, {{}, {}, {}, {}}}, {{{}, {}, \
    {}, {}}, {{}, {}, {}, {}}, {{}, {}, {},
       3.78763}, {{}, {}, {}, {}}}, {{{},
       0.848507, {}, {}}, {{}, {}, {}, {}}, {{}, {}, {}, {}}, {{}, {}, \
    {}, {}}}}
  9. 执行 Flatten[Table[
      Outer[If[Chop[#1 - #2] === 0, Chop[#1], {}] &,
       Last /@ List @@ NRoots[poly1[z, aZeros[[i]]] == 0, z],
       Last /@ List @@ NRoots[poly2[z, aZeros[[i]]] == 0, z]], {i,
       Length[aZeros]}]]
    返回 {-0.848507, -3.78763, -1.88035 - 1.88035 I, -1.88035 +
      1.88035 I, -1.21371 - 1.21371 I, -1.21371 + 1.21371 I,
     0. - 3.78763 I, 0. + 3.78763 I, 0. - 0.848507 I, 0. + 0.848507 I,
     1.21371 - 1.21371 I, 1.21371 + 1.21371 I, 1.88035 - 1.88035 I,
     1.88035 + 1.88035 I, 3.78763, 0.848507}
  10. 执行 Table[{z ->
       Flatten[(*解集合中是{a,z}对,a取自于aZeros数组z原来就是后面这一串的计算结果*)(*Sort@
         是 Sort[...] 的等价写法*)
         Outer[If[(*which are the same?*)Chop[#1 - #2] === 0,
            Chop[#1], {}] &,(*这里为Pure Function,#1 为第一个参数*)(*the zeros of \
    poly1 with respect to z*)
          Last /@ List @@
            NRoots[poly1[z, aZeros[[i]]] == 0,
             z],(*the zeros of poly2 with respect to z*)
          Last /@ List @@
            NRoots[poly2[z, aZeros[[i]]] == 0, z]]][[1]]}, {i,
      Length[aZeros]}]
    返回 {{z -> -0.848507}, {z -> -3.78763}, {z -> -1.88035 -
        1.88035 I}, {z -> -1.88035 + 1.88035 I}, {z -> -1.21371 -
        1.21371 I}, {z -> -1.21371 + 1.21371 I}, {z ->
       0. - 3.78763 I}, {z -> 0. + 3.78763 I}, {z ->
       0. - 0.848507 I}, {z -> 0. + 0.848507 I}, {z ->
       1.21371 - 1.21371 I}, {z -> 1.21371 + 1.21371 I}, {z ->
       1.88035 - 1.88035 I}, {z -> 1.88035 + 1.88035 I}, {z ->
       3.78763}, {z -> 0.848507}}
4)  将 Chop[{poly1[z, a], poly2[z, a]} /. %6] 带入两个方程中验证返回
{{0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 
  0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}}

5)数值解法求出两个方程组联立的解,执行 Sort @ Chop[NSolve[{poly1[z, a] == 0, poly2[z, a] == 0}, {a, z}]] 返回的结果与第3)节中的完全一致,由此可以验证Resultant求解方程组的正确性

{{a -> -0.697333, z -> -0.848507}, {a -> -0.682883, 

  z -> -3.78763}, {a -> -0.558864 - 0.558864 I, 

  z -> -1.88035 - 1.88035 I}, {a -> -0.558864 + 0.558864 I, 

  z -> -1.88035 + 1.88035 I}, {a -> -0.0874929 - 0.0874929 I, 

  z -> -1.21371 - 1.21371 I}, {a -> -0.0874929 + 0.0874929 I, 

  z -> -1.21371 + 1.21371 I}, {a -> 0. - 0.682883 I, 

  z -> 0. - 3.78763 I}, {a -> 0. + 0.682883 I, 

  z -> 0. + 3.78763 I}, {a -> 0. - 0.697333 I, 

  z -> 0. - 0.848507 I}, {a -> 0. + 0.697333 I, 

  z -> 0. + 0.848507 I}, {a -> 0.0874929 - 0.0874929 I, 

  z -> 1.21371 - 1.21371 I}, {a -> 0.0874929 + 0.0874929 I, 

  z -> 1.21371 + 1.21371 I}, {a -> 0.558864 - 0.558864 I, 

  z -> 1.88035 - 1.88035 I}, {a -> 0.558864 + 0.558864 I, 

  z -> 1.88035 + 1.88035 I}, {a -> 0.682883, 

  z -> 3.78763}, {a -> 0.697333, z -> 0.848507}}  


  评论这张
 
阅读(516)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2016