用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}

正确了终于。
---------------------