1. 程式人生 > >用Python解線性方程組——Scipy包和自己寫

用Python解線性方程組——Scipy包和自己寫

## 一、基於SymPy庫 
用Python解決方程組、微積分等問題,主要是用到Python的一個庫——SymPy庫。可以說這個專案也主要是學習SymPy庫的用法。

解二元一次方程功能實現

解方程的功能主要是使用Sympy中solve函式實現。

示例題目是: 


方程表示

程式碼表示與手寫還是有區別的,下面列出常用的: 
加號 + 
減號 - 
除號 / 
乘號 * 
指數 ** 
對數 log() 
e的指數次冪 exp()

題目中表達式可表示為:

2 * x - y - 3 = 0 
3 * x + y - 7 = 0

完整程式碼為:

from sympy import *
x = Symbol('x')
y = Symbol('y')
print solve([2 * x - y - 3, 3 * x + y - 7],[x, y])
 
麻麻,我跟正確答案一樣哦~ 


二、自己試圖寫個簡單的函式
求解方程組需要先消元,我隨手寫了一個這樣的方程組: 
x + 2y +3z = 6 
3x + y = 3 
-x + y = 3

2.1 先消元
按照最原始的方法消元(初中學的方法),程式碼如下:

# 寫一個n*n非奇異線性方程組求解函式
# Elimination:消元
def Elimination1(A = A):
    for i in range(1, A.shape[0]):
        A[i] = A[i] - A[0]*(A[i,0]/A[0,0])      
    return(A)

def Elimination(A = A):
    j = A.shape[0] - 1
    i = 0
    while j > 0:
        A[i:,i:] = Elimination1(A = A[i:,i:])
        i = i + 1
        j = j - 1
    return(A)
 
第一個函式是消一次,第二個函式通過呼叫第一個函式,迴圈消元。

把方程組寫進來,執行函式輸出。

import numpy as np
A = np.matrix([[1,2,3,6],[3,1,0,3],[-1,1,0,3]])
Elimination(A=A)
 
輸出結果

matrix([[ 1, 2, 3, 6], 
[ 0, -5, -9, -15], 
[ 0, 0, -2, 0]])

消元后,下一步就可以求解了。

2.2 求解方程組
程式碼有問題,求解有錯誤 
但是還沒查出來是哪裡出了問題

# 求解方程組
def SolveX(A = None):
    A = Elimination(A=A)
    X = {}
    m = A.shape[0]
    n = A.shape[1]
    for i in range(m):
        j = m-1-i
        Y = A[j,-1]
        k = i
        while(k>0):
            count = 0 # 問題在這!這樣內部迴圈count永遠不變。
            Y = Y - (A[2-k-count,-2-count]*X[2-count])
            count = count + 1
            k = k-1
        X[j] = Y/A[j,-2-i]
    return X
 
第三個迴圈錯了,明天一步一步的核對,看哪個地方錯了。 
分開一步一步核對是對的,然而成了函式又錯了,很奇怪哪裡出了問題。

修改如下:

# 寫一個n*n非奇異線性方程組求解函式
# Elimination:消元
def Elimination1(A = None):
    for i in range(1, A.shape[0]):
        A[i] = A[i] - A[0]*(A[i,0]/A[0,0])      
    return(A)

def Elimination(A = None):
    j = A.shape[0] - 1
    i = 0
    while j > 0:
        A[i:,i:] = Elimination1(A = A[i:,i:])
        i = i + 1
        j = j - 1
    return(A)

# 求解方程組
def SolveX(A = None):
    A = Elimination(A=A)
    X = {}
    m = A.shape[0]
    for i in range(m):
        j = m-1-i
        Y = A[j,-1]
        k = i
        count = 0
        while(k>0): 
            Y = Y - A[2-k-count,-2-count]*X[2-count]
            count = count + 1
            k = k-1
        X[j] = Y/A[j,-2-i]
    return X

import numpy as np
A = np.matrix([[1,2,3,6],[3,1,0,3],[-1,1,0,3]])
SolveX(A = A)
 
結果:

{0: 0.0, 1: 3.0, 2: -0.0}

正確了終於。
---------------------