【Python】常用数值方法的python实现

目录

• 解非线性方程
  • 方法综述
  • 问题分类
    • 求解一元方程
      • 解法一：SymPy.solve/nsolve函数求解
      • 解法二：迭代法
    • 求解多元方程组
      • 方法一：运用SymPy
      • 方法二：运用SciPy.optimize.fsolve()
• 解线性方程组
  • 插值法
    • 方法综述
    • 问题分类
      • 一元函数插值
        • B样条插值
      • 二元函数插值
        • 绘制2D图
        • 绘制3D图
  • 函数逼近（拟合）
  • 微分方程数值解法
    • 常微分方程
      • 方法一：SymPy.dsolve()
      • 方法二：scipy.integrate.odeint()
      • 方法三：单个函数四阶龙格-库塔法
      • 方法四：多个微分方程：欧拉法
    • 偏微分方程
• 参考链接

解非线性方程

方法综述

• 使用SymPy符号计算库中的solve()（解析解）/nsolve()（数值解）函数。
• 一元方程数值解可使用迭代法辅助求解。

问题分类

求解一元方程

问题引入：求该方程数值解

解法一：SymPy.solve/nsolve函数求解

实际上SymPy内置解方程方法为迭代法。

输入：

import numpy as np
from sympy import *
var("x")
#方法一：求解解析解，进一步得出数值解
#solve()参数：表达式，自变量
result1=solve(x**3-x-1,x)
#求出解析解后，可根据解析解找出所需要的解，evalf()输出数值解
result1[2].evalf()
#方法二：对于没有解析解的方程，直接输出数值解
result2=nsolve(ln((513+0.6651*x)/(513-0.6651*x))-x/(1400*0.0918),0)
result2

输出：

1.32471795724475
0

解法二：迭代法

本方法仅提供一种思路，实际求解时可直接用解法一。实际上SymPy内置解方程方法为迭代法。

示例：双点弦截法

输入：

#解方程法二：迭代法
#方法一：弦截法
import numpy as np
import matplotlib.pyplot as plt
#可以显示中文
plt.rcParams["font.sans-serif"] = ["SimHei"]
plt.rcParams['axes.unicode_minus'] = False
# 设置风格
plt.style.use('ggplot')
# 定义函数
init_fun = lambda x: (513+0.6651*x)/(513-0.6651*x)-x/(1400*0.0918)
# 导数
#deri_fun = lambda x: 2*x-4
# input
'''
x0:初始值1
x1:初始值2
theta:误差上界
'''
x0=float(input('输入初始点x0：较大值
'))
x1=float(input('输入初始点x1：较小值
'))
theta=1e-7
fig_1 = plt.figure(figsize = (8, 6))
plt.hlines(0,-1,x0,'black','--')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('函数图像')
# 函数图像
x=[]
if x0>0:
    x = np.arange(-1,x0,0.05)
    plt.hlines(0,-1,x0,'black','--')
else:
    x = np.arange(x0,10,0.05)
    plt.hlines(0,x0,10,'black','--')
y = init_fun(x)
def Secant(func = init_fun, x0 =x0, x1 = x1, theta = theta):
    number=0
    while True:
        x2=x0-func(x0)*(x1-x0)/(func(x1)-func(x0))
        plt.vlines(x0,0,init_fun(x0),'blue','--')
        plt.plot([x2,x0],[0,func(x0)],'r--',c='green')
        plt.scatter(x0,func(x0),c='black')
        if abs(func(x2))<theta:
            return x2,number
        x0=x1
        x1=x2
        number += 1
        
# 迭代法计算求解x0
xi,number = Secant(init_fun, x0, x1, theta)
print('迭代结果：'+str(xi))
print('迭代次数：'+str(number))
## 函数求解
plt.plot(x,y)
plt.show()

输出：

输入初始点x0：较大值
300
输入初始点x1：较小值
100
迭代结果：256.84798935747875
迭代次数：5

如果希望在lambda表达式输入对数、指数操作，只需引用np：

# 定义函数
init_fun = lambda x: np.log((513+0.6651*x)/(513-0.6651*x))-x/(1400*0.0918)

求解多元方程组

与求解一元方程同理。

方法一：运用SymPy

输入：

import numpy as np
from sympy import *
var("x,y")
#方法一：求解解析解，进一步得出数值解
#solve()参数：表达式，变量
result1=solve((x**2+x*y+1,y**2+x*y+2),x,y)
result1
#求出解析解后，可根据解析解找出所需要的解，evalf()输出数值解
result1[0][0].evalf()
#方法二：对于没有解析解的方程，直接输出数值解
var("x,y")
f1=x**2+y*x-1
f2=y**2+x+2
#初值要选好，否则会迭代不收敛
nsolve((f1,f2),(x,y),(I,I))

输出：

[(-sqrt(3)*I/3, -2*sqrt(3)*I/3), (sqrt(3)*I/3, 2*sqrt(3)*I/3)]
-0.577350269189626*I
Matrix([
[0.518912794385156 - 0.666609844932019*I],
[ 0.208223290106041 + 1.60070913439255*I]])

方法二：运用SciPy.optimize.fsolve()

输入：

from scipy.optimize import fsolve
from math import sin
def f(x):
    x0,x1,x2=x.tolist()
    return[
        5*x1+3,
        4*x0*x0-2*sin(x1*x2),
        x1*x2-1.5
    ]
#f计算方程组误差，[1,1,1]是初值
result=fsolve(f,[1,1,1])
print(result)
print(f(result))

输出：

[-0.70622057 -0.6        -2.5       ]
[0.0, -9.126033262418787e-14, 5.329070518200751e-15]

解线性方程组

方法：使用scipy.linalg的solve()方法

输入：

import numpy as np
from scipy.linalg import solve
a = np.array([[3, 1, -2], [1, -1, 4], [2, 0, 3]])
b = np.array([5, -2, 2.5])
x = solve(a, b)
x

输出：

array([0.5, 4.5, 0.5])

插值法

方法综述

插值是一种通过已知的离散数据来求未知数据的方法。与拟合不同的是，它要求曲线通过所有的已知数据。SciPy的interpolate模块提供了许多进行插值运算的函数。

问题分类

一元函数插值

其他插值方法的Python实现可参考：http://liao.cpython.org/scipy11/

文档仅演示最常用的B样条插值。

B样条插值

方法一：interp1d()

一维数据的B样条插值运算可以通过interp1d()完成。

interp1d(x,y,kind='cubic')

参数x和y：一系列已知的数据点。
kind：插值类型。

• 'zero','nearest'：阶梯插值，相当于0阶B样条曲线。
• 'slinear','linear'：线性插值，相当于1阶B样条曲线。
• 'quadratic','cubic'：2阶和3阶B样条曲线，更高阶的曲线可直接使用整数值指定。

实际常用三次B样条插值。

输入：（数据来源于本人物理实验BII弗兰克-赫兹实验数据）

#绘图工具：pylab
import numpy as np
from scipy import interpolate
import pylab as pl
x=[15.1,17.1,19.1, 19.4,21.4,23.4, 25.6,27.6,29.6, 30.6,32.6,34.6, 37.0,39.0,41.0, 42.5,44.5,46.5, 49.2,51.2,53.2, 54.7,56.7,58.7, 61.8,63.8,65.8, 67.6,69.6,71.6, 75.2,77.2,79.2, 81.2,83.2,85.2]
y=[175,213,179, 176,87,178, 298,340,272, 171,43,168, 345,407,332, 159,24,158, 399,466,371, 210,57,164, 433,517,437, 250,144,229, 518,588,505, 369,307,379]
xnew=np.arange(15.1,85.2,0.1)
pl.plot(x,y,"ro")
for kind in ["cubic"]:#插值方式
    #"nearest","zero"为阶梯插值
    #slinear 线性插值
    #"quadratic","cubic" 为2阶、3阶B样条曲线插值
    f=interpolate.interp1d(x,y,kind=kind)
    # ‘slinear’, ‘quadratic’ and ‘cubic’ refer to a spline interpolation of first, second or third order)
    ynew=f(xnew)
    pl.plot(xnew,ynew,label=str(kind))
pl.legend(loc="lower right")
pl.show()

输出：

与示波器输出的图像十分类似。

方法二：splrep()+splev()

该函数可以找到一维曲线的B-spline表示。

输入：

import numpy as np
import matplotlib.pyplot as plt
#进行样条插值
import scipy.interpolate as spi
 
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
#数据准备
X=[15.1,17.1,19.1, 19.4,21.4,23.4, 25.6,27.6,29.6, 30.6,32.6,34.6, 37.0,39.0,41.0, 42.5,44.5,46.5, 49.2,51.2,53.2, 54.7,56.7,58.7, 61.8,63.8,65.8, 67.6,69.6,71.6, 75.2,77.2,79.2, 81.2,83.2,85.2]
Y=[175,213,179, 176,87,178, 298,340,272, 171,43,168, 345,407,332, 159,24,158, 399,466,371, 210,57,164, 433,517,437, 250,144,229, 518,588,505, 369,307,379]
#定义插值点
ix3=np.arange(15.1,85.2,0.1)
 
#进行三次样条拟合
ipo3=spi.splrep(X,Y,k=3) #样本点导入，生成参数
iy3=spi.splev(new_x,ipo3) #根据观测点和样条参数，生成插值
#作图
plt.plot(X,Y,'o',ix3,iy3)
plt.xlabel('Vp/V')
plt.ylabel('IA/μA')
plt.title('弗兰克-赫兹图')
plt.show()

输出：

二元函数插值

绘制2D图

输入：

# -*- coding: utf-8 -*-
import numpy as np
from scipy import interpolate
import pylab as pl
import matplotlib as mpl
def func(x, y):
    return (x+y)*np.exp(-5.0*(x**2 + y**2))
# X-Y轴分为15*15的网格
y,x= np.mgrid[-1:1:15j, -1:1:15j]
fvals = func(x,y) # 计算每个网格点上的函数值  15*15的值
print(len(fvals[0]))
#三次样条二维插值
newfunc = interpolate.interp2d(x, y, fvals, kind='cubic')
# 计算100*100的网格上的插值
xnew = np.linspace(-1,1,100)#x
ynew = np.linspace(-1,1,100)#y
fnew = newfunc(xnew, ynew)#仅仅是y值   100*100的值
# 输出解得所求点的插值
print(newfunc(0.01,0.01))
# 绘图
# 为了更明显地比较插值前后的区别，使用关键字参数interpolation='nearest'
# 关闭imshow()内置的插值运算。
pl.subplot(121)
im1=pl.imshow(fvals, extent=[-1,1,-1,1], cmap=mpl.cm.hot, interpolation='nearest', origin="lower")#pl.cm.jet
#extent=[-1,1,-1,1]为x,y范围  favals为
pl.colorbar(im1)
pl.subplot(122)
im2=pl.imshow(fnew, extent=[-1,1,-1,1], cmap=mpl.cm.hot, interpolation='nearest', origin="lower")
pl.colorbar(im2)
pl.show()

输出：

15
[0.01985751]

绘制3D图

输入：

# -*- coding: utf-8 -*-
"""
演示二维插值。
"""
# -*- coding: utf-8 -*-
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib as mpl
from scipy import interpolate
import matplotlib.cm as cm
import matplotlib.pyplot as plt
def func(x, y):
    return (x+y)*np.exp(-5.0*(x**2 + y**2))
# X-Y轴分为20*20的网格
x = np.linspace(-1,1,20)
y = np.linspace(-1,1,20)
x, y = np.meshgrid(x, y)#20*20的网格数据
fvals = func(x,y) # 计算每个网格点上的函数值  15*15的值
fig = plt.figure(figsize=(9, 6))
#Draw sub-graph1
ax=plt.subplot(1, 2, 1,projection = '3d')
surf = ax.plot_surface(x, y, fvals, rstride=2, cstride=2, cmap=cm.coolwarm,linewidth=0.5, antialiased=True)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x, y)')
plt.colorbar(surf, shrink=0.5, aspect=5)#标注
#二维插值
newfunc = interpolate.interp2d(x, y, fvals, kind='cubic')#newfunc为一个函数
# 计算100*100的网格上的插值
xnew = np.linspace(-1,1,100)#x
ynew = np.linspace(-1,1,100)#y
fnew = newfunc(xnew, ynew)#仅仅是y值   100*100的值  np.shape(fnew) is 100*100
xnew, ynew = np.meshgrid(xnew, ynew)
# 输出解得所求点的插值
print(newfunc(0.01,0.01))
# 绘制3D图
ax2=plt.subplot(1, 2, 2,projection = '3d')
surf2 = ax2.plot_surface(xnew, ynew, fnew, rstride=2, cstride=2, cmap=cm.coolwarm,linewidth=0.5, antialiased=True)
ax2.set_xlabel('xnew')
ax2.set_ylabel('ynew')
ax2.set_zlabel('fnew(x, y)')
plt.colorbar(surf2, shrink=0.5, aspect=5)#标注
plt.show()

输出：

[0.01983083]

函数逼近（拟合）

最常用方法：最小二乘拟合

输入：

import numpy as np
from scipy.optimize import leastsq
x=np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78])
y=np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05])
def residuals(p):
    "计算以p为参数的直线和原始数据误差"
    k,b=p#若不是线性拟合，则修改对应参数
    return y-(k*x+b)
#leastsq使residuals()输出数组的平方和最小，初值[1,0]
r=leastsq(residuals,[1,0])#若不是线性拟合，则修改对应参数
k,b=r[0]#若不是线性拟合，则修改对应参数
x1=np.arange(2.66,8.71,0.01)
y1=k*x+b#若不是线性拟合，则修改对应参数
print("k=",k,"b=",b)#若不是线性拟合，则修改对应参数
#画图
import matplotlib.pyplot as plt
plt.scatter(x,y,c='g',label='scatter')#散点图
plt.plot(x1,y1,'b--',label='fitting')
plt.title('polyfitting')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()#显示标签
plt.show()

输出：

k= 0.6134953491930442 b= 1.794092543259387

计算误差曲面函数：

def S(k,b):
  #计算直线y=k*x+b与原始数据x,y误差的平方和
  error=np.zeros(k,shape)
  for  x,y in zip(x,y)
    error+=(y-(k*x+b))**2
   return error

若使用其他函数进行拟合，只需将k*x+b替换为对应和函数（参数）即可。
如，三次多项式拟合图像如下：

微分方程数值解法

常微分方程

方法一：SymPy.dsolve()

输入：

import numpy as np
from sympy import *
f = Function('f')
x = symbols('x')
eq = Eq(f(x).diff(x, x) - 2*f(x).diff(x) + f(x), sin(x))
print(dsolve(eq, f(x)))

输出：

Eq(f(x), (C1 + C2*x)*exp(x) + cos(x)/2)

方法二：scipy.integrate.odeint()

这个函数，要求微分方程必须化为标准形式，即

输入：

import math
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def func(y, t):
    return t * math.sqrt(y)
YS=odeint(func,y0=1,t=np.arange(0,10.1,0.1))
t=np.arange(0,10.1,0.1)
plt.plot(t, YS, label='odeint')
plt.legend()
plt.show()

输出：

方法三：单个函数四阶龙格-库塔法

输入：

import math
import numpy as np
import matplotlib.pyplot as plt
def runge_kutta(y, x, dx, f):
    """ y is the initial value for y
        x is the initial value for x
        dx is the time step in x
        f is derivative of function y(t)
    """
    k1 = dx * f(y, x)
    k2 = dx * f(y + 0.5 * k1, x + 0.5 * dx)
    k3 = dx * f(y + 0.5 * k2, x + 0.5 * dx)
    k4 = dx * f(y + k3, x + dx)
    return y + (k1 + 2 * k2 + 2 * k3 + k4) / 6.
if __name__=='__main__':
    t = 0.
    y = 1.
    dt = .1
    ys, ts = [], []

def func(y, t):
    return t * math.sqrt(y)

while t <= 10:
    y = runge_kutta(y, t, dt, func)
    t += dt
    ys.append(y)
    ts.append(t)
exact = [(t ** 2 + 4) ** 2 / 16. for t in ts]
plt.plot(ts, ys, label='runge_kutta')
plt.plot(ts, exact, label='exact')
plt.legend()
#plt.show()

输出：

方法四：多个微分方程：欧拉法

输入：

import numpy as np
"""
移动方程：
t时刻的位置P(x,y,z)
steps：dt的大小
sets：相关参数
"""
def move(P, steps, sets):
    x, y, z = P
    sgima, rho, beta = sets
    # 各方向的速度近似
    dx = sgima * (y - x)
    dy = x * (rho - z) - y
    dz = x * y - beta * z
    return [x + dx * steps, y + dy * steps, z + dz * steps]
# 设置sets参数
sets = [10., 28., 3.]
t = np.arange(0, 30, 0.01)
# 位置1：
P0 = [0., 1., 0.]
P = P0
d = []
for v in t:
    P = move(P, 0.01, sets)
    d.append(P)
dnp = np.array(d)
# 位置2：
P02 = [0., 1.01, 0.]
P = P02
d = []
for v in t:
    P = move(P, 0.01, sets)
    d.append(P)
dnp2 = np.array(d)
"""
画图
"""
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(dnp[:, 0], dnp[:, 1], dnp[:, 2])
ax.plot(dnp2[:, 0], dnp2[:, 1], dnp2[:, 2])
plt.show()

输出：

偏微分方程

与解常微分方程原理相同。

参考链接

• 《Python科学计算》-张若愚
• SymPy.solve()官方文档：https://docs.sympy.org/latest/tutorial/solvers.html
• SciPy.interpolate.splrep()官方文档（英文版）：https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.interpolate.splrep.html
• SciPy.interpolate.splrep()官方文档（中文版）：https://vimsky.com/examples/usage/python-scipy.interpolate.splrep.html
• Python.SciPy实现Hermite插值：http://liao.cpython.org/scipy13/
• 常微分方程数值解：Python求解：https://www.jianshu.com/p/8d3671f9148d
• 微分方程在物理学中的应用：https://zhuanlan.zhihu.com/p/81488678
  https://zhuanlan.zhihu.com/p/164627678