数值计算笔记之插值(三) 分段线性插值

2023-11-08

0、回顾:对于 拉格朗日插值多项式与牛顿插值多项式的统一

  • n次拉格朗日插值多项式为: L_{n}(x) = \sum_{i = 0}^{n}l_{i}(x)\cdot y_{i} ,其中 l_{i}(x) = \prod_{j=0,j\neq i}^{n}\frac{x-x_{j}}{x_{i}-x_{j}}
  • 牛顿插值公式:N_{n}(x) = a_{0}+a_{1}(x-x_{0})+a_{2}(x-x_{0})(x-x_{1})+\cdots +a_{n}(x-x_{0})(x-x_{1})\cdots (x-x_{n-1})

在插值节点、插值条件相同的情况下,二者本质一样,只是计算过程不一样。

牛顿插值适合需要增加节点,提高精度的情况,不需要重新开始计算,可以利用上一步的计算结果。而可以利用上一步的计算结果的特性也是牛顿插值比拉格朗日插值较好的一点。即牛顿插值具有承继性。

问题:是不是节点越多,结果精度越高呢? 肯定不是这样的。

举个例子:函数f(x) = \frac{1}{1+25x^{2}},x\in [-1,1]

1、将[-1,1] 5等分,以分点做插值节点。

x -1 -0.5 0 0.5 1
f(x) 0.0385 0.1379 1 0.1379 0.0385

为了图像化,所以使用matlab来实现牛顿法插值。代码附在最后。

2、将[-1,1] 10等分,以分点做插值节点。

x -1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1
f(x) 0.0385 0.0664 0.1379 0.39 1 0.39 0.1379 0.0664 0.0385

可以看出,当n增大时,图像在[-1,1]附近有着激烈振荡现象,这种现象称为龙格现象

matlab代码:

调用函数

function[y,A,C,L] = newtonInter(X,Y,x)
n = length(X);
m = length(x);
for t = 1 : m
    z = x(t);
    A = zeros(n,n);
    A(:,1) = Y';
        for j = 2 : n
            for i = j : n
                A(i,j) = (A(i,j-1) - A(i-1,j-1))/(X(i)-X(i-j+1));  %差商矩阵
            end
        end
        
        C = A(n, n); 
        for k = (n-1):-1:1
            C = conv(C, poly(X(k)));
            d = length(C);
            C(d) = C(d) + A(k,k);
        end
        y(t) = polyval(C,z);
end
L = poly2sym(C);

 执行文件

clear,clc;
% X = [-1, -0.5, 0, 0.5, 1];
% Y = [0.0385, 0.1379, 1, 0.1379, 0.0385];
X = [-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1];
Y = [0.0385, 0.0664, 0.1379, 0.39, 1, 0.39, 0.1379, 0.0664, 0.0385];
x = linspace(-1,1,250);  
[y,A,C,L] = newtonInter(X, Y, x);   
y1 = 1 ./ (1 + 25 * x.^2);
hold on  
plot(X, Y, 'or', x, y, '.k',x,y1,'-b');  
legend('样本点','牛顿插值估算','1/(1+25x^2)');

1、分段线性插值

对于龙格现象,计算是在整个区间上的,如果想解决这个问题,那么就把整个区间分成一段一段的来计算,这就是分段线性插值。

设 n+1 个节点

x x_{0} x_{1} \cdots x_{n}
f(x) y_{0} y_{1} \cdots y_{n}

在区间[x_{i-1},x_{i}] (i=1,2,\cdots ,n)上做线性插值,,在图形上即为把两点用线段相连,n条线段组成折线,该折线对应的函数称为分段线性插值函数。

插值基函数为:L(x) = y_{k} + (x-x_{k})\frac{y_{k+1}-y_{k}}{x_{k+1}-x_{k}}

首先要确定间隔符号k,使得x_{k}\leq x< x_{k+1}

令 s = x-x_{k}  、 \delta _{k} = \frac{y_{k+1}-y_{k}}{x_{k+1}-x_{k}} ,则 L(x) = y_{k} + (x-x_{k})\frac{y_{k+1}-y_{k}}{x_{k+1}-x_{k}} = y_{k}+s\delta _{k}

插值结果为:(节点越多,越准确)

可以看出,分段线性插值也是存在缺点的,那就是连线是折线,整体的曲线不够光滑。

matlab代码:

function v = piecelin(x,y,u)
delta = diff(y) ./ diff(x);

n = length(x);
k = ones(size(u));
for j = 2:n-1
    k( x(j) <= u) = j;
end
s = u - x(k);
v = y(k) + s.*delta(k);
end

执行代码;

clear,clc;
X = [-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1];
Y = [0.0385, 0.0664, 0.1379, 0.39, 1, 0.39, 0.1379, 0.0664, 0.0385];
x = linspace(-1,1,250);  
y = piecelin(X,Y,x);
plot(X, Y, 'or', x, y, '.k');  
legend('样本点','分段线性插值估算');

C++代码 (其计算出的点与matlab一致):

#include<iostream>
#include<vector>
#include<fstream>
using namespace std;

bool piecelin(double X[], double Y[], vector<double> &x, vector<double> &v);

int main()
{
	double x[] = { -1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1 };
	double y[] = { 0.0385, 0.0664, 0.1379, 0.39, 1, 0.39, 0.1379, 0.0664, 0.0385 };
	vector<double> u;
	for (double i = -1; i <= 1; i = i + 0.01) {
		u.push_back(i);
	}
	vector<double> v;
	/*piecelin(x, y, u,v);*/
	int  n = sizeof(x) / sizeof(x[0]);
	int k = u.size();
	for (int i = 0; i < n - 1; i++) {
		for (int j = 0; j < k; j++) {
			if (x[i] <= u[j] && u[j] < x[i + 1]) {
				double s = u[j] - x[i];
				double delta = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
				v.push_back(y[i] + s * delta);
			}

			if (u[j] > x[i + 1])
				break;
		}
	}

	int size = v.size();

	std::ofstream output;
	output.open("piecelin.txt");
	for (unsigned i = 0; i < size; i++) {
		output << u[i] << '\t' << v[i] << std::endl;
	}
	output.close();
	cout << "Hello World !" << endl;
}

 

本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

数值计算笔记之插值(三) 分段线性插值 的相关文章

  • 如何在MVVM中管理多个窗口

    我知道有几个与此类似的问题 但我还没有找到明确的答案 我正在尝试深入研究 MVVM 并尽可能保持纯粹 但不确定如何在坚持模式的同时启动 关闭窗口 我最初的想法是向 ViewModel 发送数据绑定命令 触发代码来启动一个新视图 然后通过 X
  • C# 和 Javascript SHA256 哈希的代码示例

    我有一个在服务器端运行的 C 算法 它对 Base64 编码的字符串进行哈希处理 byte salt Convert FromBase64String serverSalt Step 1 SHA256Managed sha256 new S
  • 如何在列表框项目之间画一条线

    我希望能够用水平线分隔列表框中的每个项目 这只是我用于绘制项目的一些代码 private void symptomsList DrawItem object sender System Windows Forms DrawItemEvent
  • Newtonsoft JSON PreserveReferences处理自定义等于用法

    我目前在使用 Newtonsoft Json 时遇到一些问题 我想要的很简单 将要序列化的对象与所有属性和子属性进行比较以确保相等 我现在尝试创建自己的 EqualityComparer 但它仅与父对象的属性进行比较 另外 我尝试编写自己的
  • 在 Visual Studio 2008 上设置预调试事件

    我想在 Visual Studio 中开始调试程序之前运行一个任务 我每次调试程序时都需要运行此任务 因此构建后事件还不够好 我查看了设置的 调试 选项卡 但没有这样的选项 有什么办法可以做到这一点吗 你唯一可以尝试的 IMO 就是尝试Co
  • C - 找到极限之间的所有友好数字

    首先是定义 一对友好的数字由两个不同的整数组成 其中 第一个整数的除数之和等于第二个整数 并且 第二个整数的除数之和等于第一个整数 完美数是等于其自身约数之和的数 我想做的是制作一个程序 询问用户一个下限和一个上限 然后向他 她提供这两个限
  • C 预处理器库

    我的任务是开发源分析工具C程序 并且我需要在分析本身之前预处理代码 我想知道什么是最好的图书馆 我需要一些重量轻 便于携带的东西 与其推出自己的 为什么不使用cpp这是的一部分gcc suite http gcc gnu org onlin
  • Json.NET - 反序列化接口属性引发错误“类型是接口或抽象类,无法实例化”

    我有一个类 其属性是接口 public class Foo public int Number get set public ISomething Thing get set 尝试反序列化Foo使用 Json NET 的类给我一条错误消息
  • 使用 System.Text.Json 即时格式化 JSON 流

    我有一个未缩进的 Json 字符串 例如 hash 123 id 456 我想缩进字符串并将其序列化为 JSON 文件 天真地 我可以使用缩进字符串Newtonsoft如下 using Newtonsoft Json Linq JToken
  • 如何返回 json 结果并将 unicode 字符转义为 \u1234

    我正在实现一个返回 json 结果的方法 例如 public JsonResult MethodName Guid key var result ApiHelper GetData key Data is stored in db as v
  • 在 ASP.NET Core 3.1 中使用包含“System.Web.HttpContext”的旧项目

    我们有一些用 Net Framework编写的遗留项目 应该由由ASP NET Core3 1编写的API项目使用 问题是这些遗留项目正在使用 System Web HttpContext 您知道它不再存在于 net core 中 现在我们
  • vector 超出范围后不清除内存

    我遇到了以下问题 我不确定我是否错了或者它是一个非常奇怪的错误 我填充了一个巨大的字符串数组 并希望在某个点将其清除 这是一个最小的例子 include
  • clang 实例化后静态成员初始化

    这样的代码可以用 GCC 编译 但 clang 3 5 失败 include
  • 如何使我的表单标题栏遵循 Windows 深色主题?

    我已经下载了Windows 10更新包括黑暗主题 文件资源管理器等都是深色主题 但是当我创建自己的 C 表单应用程序时 标题栏是亮白色的 如何使我自己的桌面应用程序遵循我在 Windows 中设置的深色主题 你需要调用DwmSetWindo
  • 为什么 C# Math.Ceiling 向下舍入?

    我今天过得很艰难 但有些事情不太对劲 在我的 C 代码中 我有这样的内容 Math Ceiling decimal this TotalRecordCount this PageSize Where int TotalRecordCount
  • Validation.ErrorTemplate 的 Wpf 动态资源查找

    在我的 App xaml 中 我定义了一个资源Validation ErrorTemplate 这取决于动态BorderBrush资源 我打算定义独特的BorderBrush在我拥有的每个窗口以及窗口内的不同块内
  • C 中的异或运算符

    在进行按位操作时 我在确定何时使用 XOR 运算符时遇到一些困难 按位与和或非常简单 当您想要屏蔽位时 请使用按位 AND 常见用例是 IP 寻址和子网掩码 当您想要打开位时 请使用包含或 然而 XOR 总是让我明白 我觉得如果在面试中被问
  • 如何在 C++ BOOST 中像图形一样加载 TIFF 图像

    我想要加载一个 tiff 图像 带有带有浮点值的像素的 GEOTIFF 例如 boost C 中的图形 我是 C 的新手 我的目标是使用从源 A 到目标 B 的双向 Dijkstra 来获得更高的性能 Boost GIL load tiif
  • 防止索引超出范围错误

    我想编写对某些条件的检查 而不必使用 try catch 并且我想避免出现 Index Out of Range 错误的可能性 if array Element 0 Object Length gt 0 array Element 1 Ob
  • 使用 libcurl 检查 SFTP 站点上是否存在文件

    我使用 C 和 libcurl 进行 SFTP FTPS 传输 在上传文件之前 我需要检查文件是否存在而不实际下载它 如果该文件不存在 我会遇到以下问题 set up curlhandle for the public private ke

随机推荐

  • Redis的RDB和AOF两种存储方式

    Redis其实就是一个用C语言写的一个程序 这个程序用来存储 key value数据 数据先放在内存 然后写入磁盘指定位置 这么理解十分肤浅 但tm好像就是这样啊 下面我们梳理一下Redis存储两种方式 RDB和AOF 第一种方式 RDB
  • Python求水仙花数

    水仙花数 Narcissistic Number 是一个三位数 其各位数字的立方和等于该数本身 例如 153 1 3 5 3 3 3 因此153是一个水仙花数 以下是一个简单的Python代码来找出所有的水仙花数 python def is
  • python数据分许基础

    1 单选题 单选题 下列关于数据和数据分析的说法正确的是 A 数据就是数据库中的表格 B 文字 声音 图像这些都是数据 C 数据分析的数据只能是结构化的 D 数据分析不可能预测未来几天的天气变化 正确答案 B 文字 声音 图像这些都是数据
  • 理解JavaScript的编译过程与运行机制

    JavaScript引擎 不是逐条解释执行javaScript代码 而是按照代码块一段段解释执行 所谓代码块就是使用
  • LeetCode 172. Factorial Trailing Zeroes

    Given an integer n return the number of trailing zeroes in n Follow up Could you write a solution that works in logarith
  • Cookie增删改查

    cookie 浏览器请求访问服务器 用请求头 首部行等信息来获取数据 服务器就会给出响应头 以及set cookie 给出的是唯一的cookie 之后就返回给浏览器 当浏览器第二次请求的时候 除了发送请求头之外还要发送cookie id 1
  • Python计算日照总辐射量

    import numpy as np from datetime import datetime import pandas as pd def get dayofyear date string param date string 字符串
  • YOLOv5训练目标检测数据集(小白)

    一 提前准备工作 1 利用labelimg软件给收集到的图片打标签 具体步骤网上都有 2 下载好yolov5 v6 1 源码 下载地址 https github com ultralytics yolov5 用pycharm打开 在项目目录
  • flutter 点击事件写法

    onTapUp onTapDown void onTapDown TapUpDetails details async 或者是 onTapDown e selectCommonWordIndex index mySetState
  • vue-入门篇

    1 目标 了解什么是VUE 2 vue基础 2 1 概述 官网 https cn vuejs org Vue js是一套构建用户界面的渐进式框架 Vue 采用自底向上增量开发的设计 Vue 的核心库只关注视图层 它不仅易于上手 还便于与第三
  • MySQL - 表索引概述

    索引概述 基本概念 日常生活中 我们经常会在电话号码簿中查阅 某人 的电话号码 按姓查询或者按字母排序查询 在字典中查阅 某个词 的读音和含义等等 以快速的找到特定记录 在这里 姓 和 字母 都可看作是索引 而按 姓 或者 字母 查询则是按
  • 进程控制一之进程创建、进程终止、进程等待

    进程创建 创建子进程使用fork函数 fork有两个返回值 pid t fork void pid t相当于int 失败 返回小于0的值 成功 0 返回给子进程 大于0 返回子进程的pid给父进程 fork失败原因 内存不足 创建PCB是需
  • 程序员编程艺术PDF

    程序员编程艺术 链接 https pan baidu com s 1XWk E2DIJwYRlXNGwB LHA 提取码 nptd
  • Java基础(24)——异常、处理异常的方式详解及示例

    Java基础 24 异常详解 版权声明 一 异常体系 1 概述 2 异常的根类 Throwable 3 错误 Error 4 Exception 二 异常的处理方式 1 默认的异常处理方式 2 try catch方式 1 基本知识 2 使用
  • boost::sort::block_indirect_sort相关的测试程序

    boost sort block indirect sort是Boost库中的一个排序算法 它在排序大型数据集时表现出色 本文将介绍如何使用Boost库中的block indirect sort算法 并提供一个相关的测试程序 首先 确保已经
  • 帆软设置参数框样式

    修改前 修改后 ps 取消掉参数面板的 常用参数组合 自定义初始化控件后点击文本框会弹出对应的显示框 each this options form name widgets function i item console info item
  • android.util.AndroidRuntimeException: You cannot combine custom titles with other title features

    在做项目的时候自定义一个TitleBar 但是 其中是用到TabHost ActivityGroup 左右滑动的时候 由于TabHost中有个默认的titleBar 而在哪个自己的主界面也有一个titlebar 两个冲突了所以会报错andr
  • 【算法竞赛宝典】语言之争

    算法竞赛宝典 语言之争 题目描述 代码展示 题目描述 代码展示 语言之争 include
  • linux下查看系统安装时间,Linux下如何查看系统启动时间和运行时间以及安装时间...

    1 uptime命令 输出 16 11 40 up 59 days 4 21 2 users load average 0 00 0 01 0 00 2 查看 proc uptime文件计算系统启动时间 cat proc uptime 输出
  • 数值计算笔记之插值(三) 分段线性插值

    0 回顾 对于 拉格朗日插值多项式与牛顿插值多项式的统一 次拉格朗日插值多项式为 其中 牛顿插值公式 在插值节点 插值条件相同的情况下 二者本质一样 只是计算过程不一样 牛顿插值适合需要增加节点 提高精度的情况 不需要重新开始计算 可以利用