python数值分析

python数值分析

上学期上数值分析课的时候被老师要求用python写代码，最后代码加上实验报告，写了一天终于给整完了。为了让大家不在这么煎熬秃顶，我就把我之前写的代码整理一下分享给大家。

python二分法解决方程:x^3±2*x-5

、、、

def solve_function(x):
    return x**3-2*x-5
def dichotomy(left, right,eps):
    mid = (left+right)/2
    count=0 # 统计迭代次数

while abs(solve_function(mid))>eps:
    mid = (left+right)/2
    if solve_function(left)*solve_function(mid)<=0:
        right=mid
    else:
        left=mid
    count=count+1
return count,mid

left=2
right=5
eps=0.0000001
count,middle=dichotomy(left, right,eps)
print("姓名:XXX"+"\n"+"学号:xxxxx")
print("迭代%d次得到的根是%f" %(count,middle))

最小二乘法实现价格预测

、、、

import matplotlib.pyplot as plt
import pandas as pd
import statsmodels.api as sm

from sklearn import datasets

# 导入数据源 sklearn

boston = datasets.load_boston()
X = boston.data
y = boston.target
print(X.shape)

# print(boston.DESCR)

# 输出所有属性


i = 0
while i <= 12:
    fig, ax = plt.subplots(figsize=(3, 3))
    ax.scatter(X[:, i], y, color='blue')
    i = i + 1
# 将每列数据与价格之间的关系用散点图输出

# plt.subplot(13,1,i+1) 子图实现方式

Z = pd.DataFrame(X, columns=boston.feature_names)
w = pd.DataFrame(y, columns=['MEDV'])
Z.drop('AGE', axis=1, inplace=True)
Z.drop('INDUS', axis=1, inplace=True)

# drop删除列需要手动添加axis=1

# 如果手动设定为True（默认为False），那么原数组直接就被替换。

Z_add1 = sm.add_constant(Z)
model = sm.OLS(w, Z_add1).fit()
# sm.OLS()为普通最小二乘回归模型， fit()用于拟合

print(model.summary())
print(model.params)

牛顿法求解方程:x^3-2*x-5

、、、

    def function(x):
    return x**3 - 2*x - 5
def daoshu(x):
    return 3*x**2-2

def Newton(x,eps):
    count=0
    while abs(function(x))>eps:
        x=x-function(x)/daoshu(x)
        count=count+1
    return count,x
x=3
eps=0.0000001
count,mid=Newton(x,eps)
print("姓名:xxx"+"\n"+"学号:xxx")
print("迭代%d次得到的根是%f" %(count,mid))

线性回归求解年利润

、、、

import math

def xianxinghuigui(x,y):
    N = float(len(x))
    Sx,Sy,Sx1,Sy1,Sxy=0,0,0,0,0
    for i in range(0,int(N)):
        Sx += x[i]
        Sy += y[i]
        Sx1 += x[i]*x[i]
        Sy1 += y[i]*y[i]
        Sxy += x[i]*y[i]
    a = (Sx*Sy/N-Sxy)/(Sx*Sx/N-Sx1)
    b = (Sy-Sx)/N
    k = abs(Sx*Sy/N-Sxy)/math.sqrt((Sx1-Sx*Sx/N)*(Sy1-Sy*Sy/N))
    return a,b,k
if __name__ =="__main__":
    X = [2011,2012,2013,2014,2015,2016]
    Y = [70,122,144,152,174,196,202]
    a,b,r = xianxinghuigui(X,Y)

    print("姓名:xxx"+"\n"+"学号:xxx")
    print("X=",X)
    print("Y=",Y)
    print("拟合结果：y=%10.5f x+%10.5f  r = %10.5f"%(a,b,r))
    print("2017年的利润为:%10.5f*2017+%10.5f")
    print(int(a*2017+b))
    print("2018年的利润为：")
    print(int(a*2018+b))

高斯消去法求解矩阵

、、、

import math
    	  
def funhang(c, r1, r2):
    for i in range(4):
        c[r1][i], c[r2][i] = c[r2][i], c[r1][i]

    return c


def funlie(c, r1, r2):
    for i in range(3):
        c[i][r1], c[i][r2] = c[i][r2], c[i][r1]
        x1[r1], x1[r2] = x1[r2], x1[r1]
    return c


def funfll(c, r1, r2):
    a = 0
    b = 0
    for i in range(r1, 3):
        for j in range(r1, 3):
            if r2 < abs(c[i][j]):
                r2 = abs(c[i][j])
                a = i
                b = j
    return a, b


def shuchu():
    print()
    for i in range(3):
        for j in range(4):
            print("{:.4f}".format(c[i][j]), " ", end=" ")
        print()
    print()


a1 = 0
a2 = 0
d1 = 0
d2 = 0
e = 0
c = [0.9428, 0.3475, -0.8468, 0.4127], \
    [0.3475, 1.8423, 0.4759, 1.7321], \
    [-0.8468, 0.4759, 1.2147, -0.8621]
##c = [2,-1,3,1],[4,2,5,4],[1,2,0,7]
x = [0, 0, 0]
x1 = [1, 2, 3]
print("学号：xxx  姓名：xxx")
print("原始矩阵为：")
shuchu()
for i in range(3):  # 行
    if i < 2:
        e = abs(c[i][i])
        d1, d2 = funfll(c, i, e)
        funlie(c, i, d1)
        funhang(c, i, d2)
        print("第", i + 1, "次""换主元后为：")
        shuchu()
        for j in range(i + 1, 3):  # 列
            a1 = c[j][i] / c[i][i]
            for k in range(i, 4):  # 完整一行
                # print(a,j,k)
                c[j][k] = c[j][k] - a1 * c[i][k]
        print("第", i + 1, "次""换主元后进行高斯消去：")
        shuchu()
print("使用列主元法后矩阵为：")
shuchu()
print("后向带入后的解为：")
for i in range(3):
    if (i == 0):
        x[2 - i] = c[2 - i][3] / c[2 - i][2 - i]
    else:
        for j in range(i):
            a2 = a2 + c[2 - i][2 - j] * x[2 - j]
        x[2 - i] = (c[2 - i][3] - a2) / c[2 - i][2 - i]
        a2 = 0
    print("X", x1[i], "=", x[2 - i])