-
2021-03-26 15:16:19
本文只是完成作业,重点是比较三种插值过程和结果的特点(略),主函数(略),不考虑太多健壮性,也不批量插值。创建项目
-
使用VS2019创建C++控制台应用
拉格朗日插值
公式
n次的拉格朗日插值多项式: L n ( x ) = ∑ k = 0 n [ y k ∏ i = 0 , i ≠ k n ( x − x i x k − x i ) ] L_n(x)=\displaystyle\sum_{k=0}^{n}[y_k\displaystyle\prod_{i=0,i≠k}^{n}(\frac{x-x_i}{x_k-x_i})] Ln(x)=k=0∑n[yki=0,i=k∏n(xk−xix−xi)]
其中,已知点有n+1个,其x值依次为 x 0 x_0 x0、 x 1 x_1 x1、… x i x_i xi(或 x k x_k xk)、… x n x_n xn,其y值同理。
当取n=1时,为线性插值;当取n=2时,为抛物插值。
思路
1.编程前,令上述公式中的n取n-1,于是公式为: L n − 1 ( x ) = ∑ k = 0 n − 1 [ y k ∏ i = 0 , i ≠ k n − 1 ( x − x i x k − x i ) ] L_{n-1}(x)=\displaystyle\sum_{k=0}^{n-1}[y_k\displaystyle\prod_{i=0,i≠k}^{n-1}(\frac{x-x_i}{x_k-x_i})] Ln−1(x)=k=0∑n−1[yki=0,i=k∏n−1(xk−xix−xi)]
2.这样编写的程序函数中,接收的形参包含四个部分:
- x值的数组
- y值得数组
- 点数n
- 插值位置x
3.返回参数:
- 插值结果y
4.数据类型选择:double
5.为了对比不同插值方法的效率和结果,不建立额外的“缓存”
代码
double lagrange(double arrX[], double arrY[], int n, double x) { int k, i; double temp; double y = 0; for (k = 0; k < n; k++) { temp = 1; for (i = 0; i < n; i++) { if (i != k) { temp *= ((x - arrX[i]) / (arrX[k] - arrX[i])); } } y += arrY[k] * temp; } return y; }
特点
优点:直观、简单、应用广泛。
缺点:当插值精度不够,增加新节点时,必须从头计算,不能利用已有结果。牛顿插值
为克服拉格朗日插值的缺点,实现灵活增加插值节点,以节省运算次数。
公式
N n ( x ) = f ( x 0 ) + ∑ i = 1 n [ ∏ k = 0 i − 1 ( x − x i ) ] f [ x 0 , x 1 , . . . , x i ] N_n(x)=f(x_0)+\displaystyle\sum_{i=1}^{n}[\displaystyle\prod_{k=0}^{i-1}({x-x_i})]f[x_0,x_1,...,x_i] Nn(x)=f(x0)+i=1∑n[k=0∏i−1(x−xi)]f[x0,x1,...,xi]
思路
- 根据已知数据的个数n,建立差商表,表内共需要存储 n ( n − 1 ) 2 \frac{n(n-1)}{2} 2n(n−1)个差商值。
- 根据差商表,将插值数据代入公式。
代码
计算差商表:
//计算差商表,n表示有n个节点,n>1,f为差商表数组 void calcChaShang(double x[], double y[], double f[], int n) { int i, j, k = 0; //一阶差商计算 for (i = 0; i < n - 1; i++) { f[k++] = (y[i] - y[i + 1]) / (x[i] - x[i + 1]); } //在一阶差商的基础上,计算到n-1阶 for (i = 1; i < n - 1; i++) { //j用来遍历第i阶的所有差商(共有n-i个) for (j = 0; j < n - i - 1; j++) { //调试时所用 /*double y1 = f[k - (n - i)]; double y2 = f[k - (n - i) + 1]; double x1 = x[j]; double x2 = x[j + i + 1];*/ //由于k自增会影响y[],所以y[]中的j可以去掉,否则错误写法: //f[k] = (f[k - (n - i) + j] - f[k - (n - i) + j + 1]) / (x[j] - x[j + i + 1]); f[k] = (f[k - (n - i)] - f[k - (n - i) + 1]) / (x[j] - x[j + i + 1]); k++; } } }
基于差商表的牛顿插值:
double newvalue(double* xArr, double* yArr, double* f, int n, double x) { //将连乘结果存储到数组中(multiplication) double *m = (double*)malloc(sizeof(double) * n); m[0] = 1.0; int i = 0; for (i = 0; i < n - 1; i++) { m[i + 1] = m[i] * (x - xArr[i]); } //最终计算 int k = 0;//通过此下标,访问差商数组 double y = yArr[0]; for (i = 1; i < n; i++) { y += m[i] * f[k]; k += n - i; } return y; }
线性分段插值
公式
在区间[a,b],划分 a = x 0 < x 1 < x 2 < . . . < x n = b a=x_0<x_1<x_2<...<x_n=b a=x0<x1<x2<...<xn=b,对于[a,b]之间的任一小区间 [ x j − 1 , x j ] [x_{j-1},x_j] [xj−1,xj],在该小区间上作线性插值:
f ( x ) ≈ L 1 ( x ) = y j − 1 x − x j x j − 1 − x j + y j x − x j − 1 x j − x j − 1 f(x)≈L_1(x)=y_{j-1}\frac{x-x_j}{x_{j-1}-x_j}+y_j\frac{x-x_{j-1}}{x_j-x_{j-1}} f(x)≈L1(x)=yj−1xj−1−xjx−xj+yjxj−xj−1x−xj−1
思路
- 已知数据的x值需要从小到大排列。
- 插值前先找到被插数据x值所在的分段,通过比较 x x x和 x j − 1 x_{j-1} xj−1的值,找到其所在的小区间 [ x j − 1 , x j ] [x_{j-1},x_j] [xj−1,xj]中的j的值。
- 找到j,代入公式。
代码
double linear(double xArr[], double yArr[], int n, double x) { int j; double y = 0; bool noFound = true; //如果是左边外插 if (x <= xArr[0]) { j = 0; } //如果是右边外插 else if (x >= xArr[n - 1]) { j = n - 2; } //内插 else { for (j = 1; j < n; j++) { //找到所在的分段 if (x <= xArr[j]) { j--; break; } } } y = yArr[j] * ((x - xArr[j + 1]) / (xArr[j] - xArr[j + 1])) + yArr[j + 1] * ((x - xArr[j]) / (xArr[j + 1] - xArr[j])); return y; }
参考
公式来自《数值计算方法及其应用》(朱长青编著,科学出版社)
更多相关内容 -
-
Python实现分段线性插值
2021-01-21 18:36:14本文实例为大家分享了Python实现分段线性插值的具体代码,供大家参考,具体内容如下 函数: 算法 这个算法不算难。甚至可以说是非常简陋。但是在代码实现上却比之前的稍微麻烦点。主要体现在分段上。 图像效果 ... -
分段线性插值Python实现(同时估计误差)
2018-05-27 16:34:13函数 y=11+x2y=11+x2y = \frac{1}{1+x^2} 算法 这个算法不算难。...主要体现在分段上。 代码 import numpy as np from sympy import * import matplotlib.pyplot as plt def f(x): return 1 / (1 ...函数
y=11+x2 y = 1 1 + x 2算法
这个算法不算难。甚至可以说是非常简陋。但是在代码实现上却比之前的稍微麻烦点。主要体现在分段上。
图像效果
代码
import numpy as np from sympy import * import matplotlib.pyplot as plt def f(x): return 1 / (1 + x ** 2) def cal(begin, end): by = f(begin) ey = f(end) I = (n - end) / (begin - end) * by + (n - begin) / (end - begin) * ey return I def calnf(x): nf = [] for i in range(len(x) - 1): nf.append(cal(x[i], x[i + 1])) return nf def calf(f, x): y = [] for i in x: y.append(f.subs(n, i)) return y def nfSub(x, nf): tempx = np.array(range(11)) - 5 dx = [] for i in range(10): labelx = [] for j in range(len(x)): if x[j] >= tempx[i] and x[j] < tempx[i + 1]: labelx.append(x[j]) elif i == 9 and x[j] >= tempx[i] and x[j] <= tempx[i + 1]: labelx.append(x[j]) dx = dx + calf(nf[i], labelx) return np.array(dx) def draw(nf): plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False x = np.linspace(-5, 5, 101) y = f(x) Ly = nfSub(x, nf) plt.plot(x, y, label='原函数') plt.plot(x, Ly, label='分段线性插值函数') plt.xlabel('x') plt.ylabel('y') plt.legend() plt.savefig('1.png') plt.show() def lossCal(nf): x = np.linspace(-5, 5, 101) y = f(x) Ly = nfSub(x, nf) Ly = np.array(Ly) temp = Ly - y temp = abs(temp) print(temp.mean()) if __name__ == '__main__': x = np.array(range(11)) - 5 y = f(x) n, m = symbols('n m') init_printing(use_unicode=True) nf = calnf(x) draw(nf) lossCal(nf)
-
[Python] 分段线性插值
2020-12-08 13:53:25利用线性函数做插值每一段的线性函数:#Program 0.6 Linear Interploationimport numpy as npimport matplotlib.pyplot as plt#分段线性插值闭包def get_line(xn, yn):def line(x):index = -1#找出x所在的区间for i ...利用线性函数做插值
每一段的线性函数:
#Program 0.6 Linear Interploation
import numpy as np
import matplotlib.pyplot as plt
#分段线性插值闭包
def get_line(xn, yn):
def line(x):
index = -1
#找出x所在的区间
for i in range(1, len(xn)):
if x <= xn[i]:
index = i-1
break
else:
i += 1
if index == -1:
return -100
#插值
result = (x-xn[index+1])*yn[index]/float((xn[index]-xn[index+1])) + (x-xn[index])*yn[index+1]/float((xn[index+1]-xn[index]))
return result
return line
xn = [i for i in range(-50,50,10)]
yn = [i**2 for i in xn]
#分段线性插值函数
lin = get_line(xn, yn)
x = [i for i in range(-50, 40)]
y = [lin(i) for i in x]
#画图
plt.plot(xn, yn, 'ro')
plt.plot(x, y, 'b-')
plt.show()
-
cpu-temp-approximation:一个使用分段线性插值和三次样条曲线(贝塞尔曲线)近似时间序列数据的类项目
2021-02-15 22:07:16_____ ______ _ _ _____ / __ \| ___ \ | | | |_ _| | / \/| |_/ / | | | | | ___ _ __ ___ _ __ | | | __/| | | | | |/ _ \ '_ ` _ \| '_ \ | \__/\| | | |_| | | | __/ | | | | | |_) | ... -
分段线性插值——数值计算方法
2011-09-24 12:34:20分段线性插值——数值计算方法 分段线性插值——数值计算方法 分段线性插值——数值计算方法 分段线性插值——数值计算方法 -
龙格现象及分段线性插值
2020-09-26 01:06:34龙格现象及分段线性插值 python画图代码 import numpy as np import matplotlib.pyplot as plt def Lagrange(arr_x, arr_y, _x): l = [0 for j in range(len(arr_x))] result = 0 for i in range(0, len(arr...龙格现象及分段线性插值
python画图代码
import numpy as np import matplotlib.pyplot as plt def Lagrange(arr_x, arr_y, _x): l = [0 for j in range(len(arr_x))] result = 0 for i in range(0, len(arr_x)): denominator = 1 molecular = 1 for j in range(0, len(arr_x)): if i != j: denominator = denominator * (arr_x[i] - arr_x[j]) molecular = molecular * (_x - arr_x[j]) l[i] = molecular / denominator result = result + l[i] * arr_y[i] return result original_x = np.arange(-5.0, 5.01, 0.01) original_y = [0.0 for j in range(len(original_x))] for i in range(len(original_y)): original_y[i] = 1 / (1 + original_x[i] * original_x[i]) x_arr = np.arange(-5.0, 5.5, 1) y_arr = [0.0 for i in range(len(x_arr))] for i in range(len(x_arr)): y_arr[i] = 1 / (1 + x_arr[i] ** 2) x = np.arange(-5.0, 5.01, 0.01) y = [0.0 for j in range(len(x))] for i in range(len(y)): y[i] = Lagrange(x_arr, y_arr, x[i]) plt.plot(original_x, original_y, label='f(x) = 1 / (1 + x2)') plt.scatter(x_arr, y_arr, label='The interpolation points') plt.plot(x, y, label='Lagrange interpolation') plt.plot([-5.5, 5.5], [0, 0], linestyle='--') plt.plot([0, 0], [-0.5, 2], linestyle='--') plt.plot(x_arr, y_arr, linestyle='--', label='piecewise linear interpolation') plt.title("Runge phenomenon, piecewise linear interpolation") plt.legend(loc="lower left") plt.xlabel("x") plt.ylabel("y") plt.show()
-
数学建模准备 插值(拉格朗日多项式插值,牛顿多项式插值,分段线性插值,分段三次样条插值,分段三次...
2019-09-03 23:20:46文章目录一、基础概念插值是什么拟合是什么插值和拟合的相同点插值和拟合的不同点二、常用的基本插值方法高次多项式插值拉格朗日多项式插值牛顿插值差商矩阵低次多项式插值(不易震荡)分段线性插值Hermite插值三次... -
计算方法(3)——分段插值法(附Python程序)
2020-12-12 13:05:39这一讲主要分段插值中的分段线性插值和分段Hermite插值,并给出分段插值的Python程序。在此之前需要注意一下,n为区间数,n+1为插值节点的个数。分段线性插值分段线性插值,需要两个列表,一个用于存放各点的x坐标,... -
numpy-如何在Python中应用分段线性拟合?
2020-12-15 12:32:46numpy-如何在Python中应用分段线性拟合?我正在尝试为数据集拟合如图1所示的分段线性拟合该图是通过设置在线获得的。 我试图使用代码来应用分段线性拟合:from scipy import optimizeimport matplotlib.pyplot as ... -
【数值分析实验】插值与拟合:拉格朗日插值、牛顿插值、分段插值;线性拟合、最小二乘拟合(python)
2021-12-29 22:01:23插值与拟合插值拉格朗日插值牛顿插值分段Hermite插值拟合最小二乘法拟合例题最小二乘拟合变式例题 调包 import math import numpy as np import sympy as sp import matplotlib.pyplot as plt 插值 拉格朗日插值 #... -
Python线性插值
2021-10-18 15:55:17Python线形插值案例一:案例二:案例三: 原文链接:https://www.jianshu.com/p/48b48e6b67cd 在算法分析过程中,我们经常会遇到数据需要处理插值的过程,为了方便理解,... -
ANSYS FLUENT 超临界流体物性分段线性插值数据批量导入
2021-12-13 09:21:37使用fluent自带数据库进行超临界流体仿真时,时不时会出现计算错误,考虑是数据库不稳定引起,而使用自定义的分段线性插值物性则不会出现上述情况。 超临界流体在拟临界点附近物性变化剧烈,需要使用较多的数据点... -
Python实现线性插值和三次样条插值
2020-11-29 09:05:13(1)、函数y = sin(x)(2)、数据准备#数据准备X=np.arange(-np.pi,np.pi,1) #定义样本点X,从-pi到pi每次间隔1Y= np.sin(X)#定义样本点Y,形成sin函数new_x=np.arange(-np.pi,np.pi,0.1) #定义差值点(3)、样条插值#... -
python线性插值解析
2021-01-14 07:46:36python线性插值解析在缺失值填补上如果用前后的均值填补中间的均值,比如,0,空,1,我们希望中间填充0.5;或者0,空,空,1,我们希望中间填充0.33,0.67这样。可以用pandas的函数进行填充,因为这个就是线性插值...