数学建模学习（1）之插值与拟合篇

小吴老师929

已于 2022-08-11 17:38:07 修改

阅读量648

点赞数

分类专栏：数学建模文章标签： python

于 2022-08-10 16:02:27 首次发布

本文链接：https://blog.csdn.net/sagewbx/article/details/126268407

版权

数学建模专栏收录该内容

1 篇文章 0 订阅

订阅专栏

目录

前言

一、插值

1.一维插值

2.二维插值

二、拟合

总结

引用

前言

注：本文仅用于自我学习，如有错误，欢迎沟通交流

下载了司老师的《python数学实验与建模》，发现比matlab版本可读性高的很多。开始了我的数学建模国赛冲刺之路！立个flag，国赛之前把这本书给刷完！冲冲冲

本章的学习要求：掌握插值和拟合的方法以及适用条件

插值与拟合的定义（参考知乎答主莫大枪）：

插值：构造已知离散点列，完全经过点列的函数

比如：三次样条插值

拟合：构造已知离散点列，整体上靠近点列的函数

比如：最小二乘法

一、插值

1.一维插值

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

实现线性插值与三次样条插值代码如下：

#一维插值
    #例见python数学实验与建模p221
x=np.arange(0,25,2)#输入x轴 include start,exclude stop,2 means step(步长）
y=np.array([12,9,9,10,18,24,28,27,25,20,18,15,13])#y值
xnew=np.linspace(0,24,500)#需修改？
f1=interp1d(x,y);y1=f1(xnew);#x，y后面不加意味着线性
f2=interp1d(x,y,'cubic');y2=f2(xnew);#cubic代表三次样条
plt.rc('font',size=16);plt.rc('font',family='SimHei')
plt.subplot(121),plt.plot(xnew,y1);plt.xlabel("(A)分段线性插值")
plt.subplot(122),plt.plot(xnew,y2);plt.xlabel("(B)三次样条插值")
plt.savefig("figure7_4.png",dpi=500);plt.show()

2.二维插值

分为二维网格节点插值和二维散乱点插值（书p222中前者使用三次样条插值，后者使用邻点插值）

参考引用【2】

x = np.arange(0, 1500, 100)
y = np.arange(1200, -100, -100)

在这里本来x是0到1400的，因为stop不包括，则额外加一个100；而且因为x要和y对应起来，y是从1200开始。

二、拟合

拟合函数的选择：

数据分布接近直线，宜选用线性函数 $f(x)=a_{1}x+a_{2}$ 拟合；

数据分布接近于抛物线，宜选用二次多项式 $f(x)=a_{1}x^{2}+a_{2}x+a_{3}$ 拟合；

如果数据分布特点是开始上升较快随后逐渐变缓，则宜选用双曲线型函数 $f(x)=\frac{x}{a_{1}x+a_{2}}$ 或指数型函数 $f(x)=a_{1}e^{-\frac{a_{2}}{x}}$ 拟合(指数型函数见下图1）；

如果数据分布特点是开始下降较快随后逐渐变缓，则宜选用 $f(x)=\frac{1}{a_{1}x+a_{2}}$ 或 $f(x)=\frac{1}{a_{1}x^{2}+a_{2}}$ 或 $f(x)=a_{1}e^{-a_{2}x}$ 等函数拟合。

常被选用的非线性拟合函数有对数函数 $f(x)=a_{1}+a_{2}lnx$ ，S型曲线函数 $f(x)=\frac{1}{a+be^{-x}}$ 等。

开始上升慢，后面上升快的用指数函数？为啥书中没有提到呢？

图一

polyfit用法

from numpy import polyfit,polyval,array,arange
from matplotlib.pyplot import plot,show,rc
x0=arange(0,1.1,0.1)
y0=array([-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2])
p=polyfit(x0,y0,2)#拟合二次多项式
print("拟合二次多项式的从高次幂到低次幂系数分别：",p)
yhat=polyval(p,[0.25,0.35]);print("预测值分别为：",yhat)
rc('font',size=16)
plot(x0,y0,'*',x0,polyval(p,x0),'-');show()

curve_fit用法

import numpy as np
from scipy.optimize import curve_fit
from matplotlib.pyplot import plot, show, rc

#调用格式为 popt,cov=curve_fit(func,xdata,ydata)
y = lambda x, a, b, c: a * x ** 2 + b * x + c#设置function
x0 = np.arange(0, 1.1, 0.1)
y0 = np.array([-0.447, 1.978, 3.28, 6.16, 7.08, 7.34, 7.66, 9.56, 9.48, 9.30, 11.2])
popt, pcov = curve_fit(y, x0, y0)
print("拟合的参数值为：", popt)
print("预测值分别为：", y(np.array([0.25, 0.35]), *popt))#预测x=0.25，0.35时候的值

rc('font', size=16)
plot(x0, y0, '*', x0, y(np.array(x0), *popt), '-')
show()

三、练习一

插值处理没有问题，问题在于三次样条插值之后如何做积分，我尝试过积分，但是样条实际上还是不连接的，只是视觉上连接，不可以有效的插值，一百个样条会有一百个积分结果，或许可以试试相加起来？

自己尝试的三次样条＋积分代码如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from scipy.integrate import quad

#做积分，哪种积分方法更好呢？
def trapezoid(f, n, a, b):  # 定义梯形公式的函数
    xi = np.linspace(a, b, n)
    h = (b - a) / (n - 1)
    return h * (np.sum(f(xi)) - (f(a) + f(b)) / 2)

def simpson(f, n, a, b):  # 定义辛普森公式的函数
    xi, h = np.linspace(a, b, 2 * n + 1), (b - a) / (2.0 * n)
    xe = [f(xi[i]) for i in range(len(xi)) if i % 2 == 0]
    xo = [f(xi[i]) for i in range(len(xi)) if i % 2 != 0]
    return h * (2 * np.sum(xe) + 4 * np.sum(xo) - f(a) - f(b)) / 3.0

x=np.array([0,2,4,5,6,7,8,9,10.5,11.5,12.5,14,16,17,18,19,20,21,22,23,24])
y=np.array([2,2,0,2,5,8,25,12,5,10,12,7,9,28,22,10,9,11,8,9,3])
xnew=np.linspace(0,24,100)#输入x轴strat and stop，以及x轴长度
f1=interp1d(x,y,'cubic');y1=f1(xnew);#cubic代表三次样条
plt.rc('font',size=16);plt.rc('font',family='simhei')
plt.subplot(121),plt.plot(xnew,y1);
plt.xlabel("三次样条插值")
plt.savefig("figureex7_1.png",dpi=500);plt.show()

a=0
b=24
n=1000
f = lambda xnew: y1

print("梯形积分I1=", trapezoid(f,n, a, b))
print("辛普森积分I2=", simpson(f, n, a, b))
print("SciPy积分I3=", quad(f, a, b))

运行结果如下：