Python求解数值分析问题（线性方程组、非线性方程、函数插值、最小二乘法、数值积分）

目录

 前言

一、线性方程组      

二、非线性方程

1、二分法

2、牛顿法 

三、函数插值

四、最小二乘法

五、数值积分

---

 前言

NumPy 和 SciPy 是 Python 进行科学计算和数据分析非常重要的两个基础包,它们各有侧重。

NumPy: 提供高性能的多维数组对象 ndarray 和丰富的数组操作方法；包含基础数学、统计、线性代数、随机数生成等功能；数组广播、向量化操作,提高计算效率；与其他Python科学计算包集成良好。

SciPy: 在NumPy基础上提供更多数学算法和高级数据处理工具；包含优化、积分、插值、拟合、特殊函数等模块；提供信号处理、图像处理、常微分方程求解等功能；统计分析模块,如分布拟合、假设检验等；提供有效的数值计算算法,可解决多种数学问题

总结:NumPy提供数组和基础运算,侧重数学计算；SciPy在NumPy基础上提供更高级的科学与工程计算算法,侧重数据处理与分析；两者DAILY配合使用,构成 Python强大的科学计算生态系统。

---

 

一、线性方程组      

 例1

 

 代码如下：

import numpy as np
# a=np.mat('1,2,3;2,4,8;9,6,3')
a=np.mat('-11,2,4,3,1;1,18,6,-3,5;3,-1,9,1,-2;3,-2,4,13,2;4,-2,3,-8,18')
b=np.mat('-9;4;-12;6;-15')
# b=np.mat('1;1;3')
c=np.linalg.solve(a,b)
print(c)

结果截图：

二、非线性方程

例2

 

 

1、二分法

代码如下：

def f(x):
    return x**5-3*x**4+6*x**3-x**2+7*x-2

# 确定区间并设置精度
left = 0
right = 1
epsilon = 0.0001

# 二分法求解
while right - left > epsilon:
    middle = (left + right) / 2  # 计算中间点
    if f(middle) == 0:
        break  # 中间点是根,结束循环
    if f(left) * f(middle) < 0: # 根处于(left, middle)区间
        right = middle
    else:   # 根处于(middle, right)区间
        left = middle

# 输出结果
print((left + right) / 2)

结果截图： 

2、牛顿法 

代码如下：

def f(x):
    return x**5-3*x**4+6*x**3-x**2+7*x-2

# 定义函数的导数
def f_prime(x):
    return 5*x**4-12*x**3+18*x**2+7

# 设置初始值和精度
x0 = 0
epsilon = 0.0001

# 牛顿法迭代
x1 = x0 - f(x0)/f_prime(x0)
while abs(x1 - x0) > epsilon:
    x0 = x1
    x1 = x0 - f(x0)/f_prime(x0)

# 输出结果
print(x1)

结果截图： 

三、函数插值

例3

x	9	11	13	15	17
lnx	2.197225	2.397895	2.564949	2.708050	2.833213

求出ln11.63 ,   ln14.13.           (2.4536)    （2.6483）

代码如下：

def f(x):
      return np.log(x)
# 构建差商表
def generate_diff_table(x_data, y_data):
    n = len(x_data)
    diff_table = [y_data]
    for r in range(1, n):
        diff_row = [0] * (n - r)
        for j in range(n - r):
            diff_row[j] = (diff_table[r-1][j+1] - diff_table[r-1][j]) / (x_data[j+r] - x_data[j])
        diff_table.append(diff_row)
    return diff_table

# 插值求解
def interpolation(x_data, x):
    n = len(x_data)
    y = 0
    diff_table = generate_diff_table(x_data, [f(x_i) for x_i in x_data])

    for i in range(n):
        p = 1
        for j in range(i):
            p = p * (x - x_data[j])
        y += diff_table[i][0] * p

    return y

x_data = [9, 11, 13, 15, 17]
x = 11.63
result = interpolation(x_data, x)
print(result)

y_data = [9, 11, 13, 15, 17]
y = 14.13
result = interpolation(y_data, y)
print(result)

结果截图： 

四、最小二乘法

例4

x	26.8	25.4	28.9	27.7	23.9	24.7	28.1	26.9	27.4	22.6	25.6
y	26.5	27.3	24.2	23.6	25.9	26.3	22.5	21.7	21.4	25.8	24.9

代码如下：

# 拟合曲线
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
from scipy.optimize import leastsq

# 样本数据
Xi = np.array([26.8, 25.4, 28.9, 27.7, 23.9, 24.7, 28.1, 26.9, 27.4, 22.6,25.6])
Yi = np.array([26.5, 27.3, 24.2, 23.6, 25.9, 26.3, 22.5, 21.7, 21.4, 25.8, 24.9])


# 需要拟合的函数func（）指定函数的形状
def func(p, x):
    k, b = p
    return k*x + b


# 定义偏差函数，x，y为数组中对应Xi,Yi的值
def error(p, x, y):
    return func(p, x) - y


# 设置k，b的初始值，可以任意设定，经过实验，发现p0的值会影响cost的值：Para[1]
p0 = [1, 20]

# 把error函数中除了p0以外的参数打包到args中,leastsq()为最小二乘法函数
Para = leastsq(error, p0, args=(Xi, Yi))
# 读取结果
k, b = Para[0]
print('k=', k, 'b=', b)

# 画样本点
plt.figure(figsize=(8, 6))
plt.scatter(Xi, Yi, color='red', label='Sample data', linewidth=2)

# 画拟合直线
x = np.linspace(30, 20, 10)
y = k * x + b

# 绘制拟合曲线
plt.plot(x, y, color='blue', label='Fitting Curve', linewidth=2)
plt.legend()  # 绘制图例

plt.xlabel('X-value', fontproperties='simHei', fontsize=12)
plt.ylabel('Y-value', fontproperties='simHei', fontsize=12)

plt.show()

结果截图：

五、数值积分

例5

代码如下： 

import numpy as np
x = np.linspace(0, 3, 1001)
# f = lambda x: x**3 - 4*x**2 + 4*x + 2
f = lambda x: 1/(1+x)
a = 0
b = 2
Ax = np.linspace(a, b, 101)
Ay = f(Ax)
def defInt_middle(f, a, b, N):
    # middle point
    result = 0; FX = []; Xn = []
    dx = abs(b - a)/N
    while a < b:
        result += f(a + dx/2)*dx
        FX += [f(a + dx/2)]
        Xn += [a]
        a += dx
    return result, FX, Xn, dx

N = 4
I_mid, FX, Xn, dx = defInt_middle(f, a, b, N)
print(I_mid)

结果截图：