2012-09-13 77 views
2

我正在嘗試使用python來解決6個非線性方程組的系統。有9個變量,其中3個是預先確定的(留下6個未知方程的6個方程組)。問題是,它可能是任何3,我無法事先知道。將預定參數傳遞給scipy.optimize.fsolve

以下是方程式(如果您有興趣)。

C11 * C12 + C21·C22 + C31·C32 = 0

C11 * C13 + C21·C23 + C31·C33 = 0

C12 * C13 + C22·C23 + C32·C33 = 0

C11 * C21 + C12·C22 + C13·C23 = 0

C11 * C31 + C12·C32 + C13·C33 = 0

C21 * C31 + C22·C32 + C23·C33 = 0


注:這是我認爲最快/最容易解決的方式。另一種可能的表達式爲:

|c11 c21 c31| 
A = |c12 c22 c32| 
    |c13 c23 c33| 

    |c11 c12 c13| 
B = |c21 c22 c23| 
    |c31 c32 c33| 

     |1 0 0| 
A*B = |0 1 0| 
     |0 0 1| 

我的問題是:反正是有設置的這些3爲固定的,並具有scipy.optimize.fsolve(或更適當的模塊?)求解剩餘參數?

回答

0

所以,我自己找到了一個有效的解決方案。不知道這是否是最好的解決方案,但它的功能。

要回答我的問題,scipy.optimize.fsolve需要一個參數args =(這裏是額外的參數)。我把預定的參數放在這裏。當函數被調用時,args首先被解析並且3個預定值被放置在適當的點上。

剩餘的6個變量被傳遞到列表中,並被迭代以填充剩餘的空白。由於參數沒有改變,每個變量總是放在矩陣中的同一個點上。

使用此方法,可以預先確定任何3個矩陣元素,並且fsolve將嘗試確定餘數。

的fsolve調用語句如下:

paramSolve1, infodict, ier, mesg = scipy.optimize.fsolve(func,(i,i,i,i,i,i),args = (knownVals[0],knownVals[1],knownVals[2]), full_output = True, warning = False) 

knwonVals是預定的參數列表,而我是一個開始猜測(全6個人失蹤參數得到了同一起跑線猜測)。 full_output允許返回可選輸出,並且當找不到解決方案時,warning = False會關閉警告消息。欲瞭解更多信息,請查看http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html

對於那些有興趣的人,問題的整個代碼如下。

import scipy 
from scipy.optimize import fsolve 

def func(params, *args): 
    c = propMatrix(createMatrix(args), params) 

    ans =(scipy.dot(c[:, 0],c[:, 1]), scipy.dot(c[:, 1],c[:, 2]), scipy.dot(c[:, 0],c[:, 2]),scipy.dot(c[:, 0],c[:, 0])-1,scipy.dot(c[:, 1],c[:, 1])-1,scipy.dot(c[:, 2],c[:, 2])-1) 
    return ans 

def createMatrix(knownVals): 

    c = [['____', '____', '____'],['____', '____', '____'], ['____', '____', '____']] 

    for element in knownVals: 
     x, y, val = element 
     c[y][x] = float(val) 
    return c 

def propMatrix(c, params): 
    for p in params: 
     assign = True 
     for x in range(3): 
      for y in range(3): 
       if c[x][y]=='____' and assign: 
        c[x][y] = float(p) 
        assign = False 

    return scipy.array(c) 

def test(c): 
    v1 = c[:, 0] 
    v2 = c[:, 1] 
    v3 = c[:, 2] 
    h1 = c[0, :] 
    h2 = c[1, :] 
    h3 = c[2, :] 
    ans = (scipy.dot(v1,v1)-1, scipy.dot(v1,v2), scipy.dot(v1, v3), scipy.dot(v2, v2)-1, scipy.dot(v2, v3), scipy.dot(v3,v3)-1, scipy.dot(h1,h1)-1, scipy.dot(h1,h2), scipy.dot(h1, h3), scipy.dot(h2, h2)-1, scipy.dot(h2, h3), scipy.dot(h3,h3)-1) 
    return ans 

def getInput(): 
    knownVals = [] 
    print """\n\nThis module analytically solves for the rotation matrix\n 
    First, enter 3 known values of the matrix:\n 
       x 
      1 2 3 
     1 | c11 c12 c13 | 
    y 2 | c21 c22 c23 | 
     3 | c31 c32 c33 |\n\n""" 

    for i in range(3): 
     invalid = True 
     print "Point Number %i:"%(i) 
     while invalid: 
      x = int(raw_input("\tx-coordinate:"))-1 
      if x>2 or x<0: 
       print "\tInvalid x-coordinate." 
      else: 
       invalid = False 
     invalid = True 
     while invalid: 
      y = int(raw_input("\ty-coordinate:"))-1 
      if y>2 or y<0: 
       print "\tInvalid y-coordinate." 
      else: 
       invalid = False 
     invalid = True 
     while invalid: 
      val = float(raw_input("\tValue:")) 
      if val>1 or val<-1: 
       print "\tInvalid value. Must be -1 <= value <= 1" 
      else: 
       invalid = False 
     knownVals.append((x, y, val)) 
    c = createMatrix(knownVals) 
    print "Input Matrix:\n\n", scipy.array(c) 
    choice = raw_input("\nIs this correct (y/n)? ") 
    if choice == "y": 
     return knownVals 
    elif choice == "n": 
     return getInput() 

def Main(): 
    solution = False 
    knownVals = getInput() 
    for i in (-1,-.5,0,.5,1): 
     paramSolve1, infodict, ier, mesg = scipy.optimize.fsolve(func,(i,i,i,i,i,i),args = (knownVals[0],knownVals[1],knownVals[2]), full_output = True, warning = False) 
     if ier == 1: 
      print "\nInitial value: %r"%(i) 
      print propMatrix(createMatrix(knownVals),paramSolve1) 
      solution = True 
    if not solution: 
     print "Could not find a valid solution" 

scipy.set_printoptions(precision = 4, suppress = True) 
Main()