详解Python牛顿插值法
作者:C-S=Cong 发布时间:2023-03-05 05:58:27
一、牛顿多项式
拉格朗日多项式的公式不具备递推性,每个多项式需要单独构造。但很多时候我们需要从若干个逼近多项式选择一个。这个时候我们就需要一个具有递推关系的方法来构造——牛顿多项式
这里的的a0,a1…等可以通过逐一带入点的值来求得。但是当项数多起来时,会发现式子变得很大,这个时候我们便要引入差商的概念(利用差分思想)具体见式子能更好理解
这里在编程实现中我们可以推出相应的差商推导方程
d(k,0)=y(k)
d(k,j)=(d(k,j-1)-d(k-1,j-1)) / (x(k)-x(k-j))
二、例题
【问题描述】考虑[0,3]内的函数y=f(x)=cos(x)。利用多个(最多为6个)节点构造牛顿插值多项式。
【输入形式】在屏幕上依次输入在区间[0,3]内的一个值x*,构造插值多项式后求其P(x*)值,和多个节点的x坐标。
【输出形式】输出牛顿插值多项式系数向量,差商矩阵,P(x*)值(保留6位有效数字),和与真实值的绝对误差(使用科学计数法,保留小数点后6位有数字)。
【样例1输入】
0.8
0 0.5 1
【样例1输出】
-0.429726
-0.0299721
1
1 0 0
0.877583 -0.244835 0
0.540302 -0.674561 -0.429726
0.700998
4.291237e-03
【样例1说明】
输入:x为0.8,3个节点为(k, cos(k)),其中k = 0, 0.5, 1。
输出:
牛顿插值多项式系数向量,表示P2(x)=-0.429726x^2 - 0.0299721x + 1;
3行3列的差商矩阵;
当x为0.8时,P2(0.8)值为0.700998
与真实值的绝对误差为:4.291237*10^(-3)
【评分标准】根据输入得到的输出准确
三、ACcode:
C++(后面还有python代码)
/*
* @Author: csc
* @Date: 2021-04-30 08:52:45
* @LastEditTime: 2021-04-30 11:57:46
* @LastEditors: Please set LastEditors
* @Description: In User Settings Edit
* @FilePath: \code_formal\course\cal\newton_quo.cpp
*/
#include <bits/stdc++.h>
#define pr printf
#define sc scanf
#define for0(i, n) for (i = 0; i < n; i++)
#define for1n(i, n) for (i = 1; i <= n; i++)
#define forab(i, a, b) for (i = a; i <= b; i++)
#define forba(i, a, b) for (i = b; i >= a; i--)
#define pb push_back
#define eb emplace_back
#define fi first
#define se second
#define int long long
#define pii pair<int, int>
#define vi vector<int>
#define vii vector<vector<int>>
#define vt3 vector<tuple<int, int, int>>
#define mem(ara, n) memset(ara, n, sizeof(ara))
#define memb(ara) memset(ara, false, sizeof(ara))
#define all(x) (x).begin(), (x).end()
#define sq(x) ((x) * (x))
#define sz(x) x.size()
const int N = 2e5 + 100;
const int mod = 1e9 + 7;
namespace often
{
inline void input(int &res)
{
char c = getchar();
res = 0;
int f = 1;
while (!isdigit(c))
{
f ^= c == '-';
c = getchar();
}
while (isdigit(c))
{
res = (res << 3) + (res << 1) + (c ^ 48);
c = getchar();
}
res = f ? res : -res;
}
inline int qpow(int a, int b)
{
int ans = 1, base = a;
while (b)
{
if (b & 1)
ans = (ans * base % mod + mod) % mod;
base = (base * base % mod + mod) % mod;
b >>= 1;
}
return ans;
}
int fact(int n)
{
int res = 1;
for (int i = 1; i <= n; i++)
res = res * 1ll * i % mod;
return res;
}
int C(int n, int k)
{
return fact(n) * 1ll * qpow(fact(k), mod - 2) % mod * 1ll * qpow(fact(n - k), mod - 2) % mod;
}
int exgcd(int a, int b, int &x, int &y)
{
if (b == 0)
{
x = 1, y = 0;
return a;
}
int res = exgcd(b, a % b, x, y);
int t = y;
y = x - (a / b) * y;
x = t;
return res;
}
int invmod(int a, int mod)
{
int x, y;
exgcd(a, mod, x, y);
x %= mod;
if (x < 0)
x += mod;
return x;
}
}
using namespace often;
using namespace std;
int n;
signed main()
{
auto polymul = [&](vector<double> &v, double er) {
v.emplace_back(0);
vector<double> _ = v;
int m = sz(v);
for (int i = 1; i < m; i++)
v[i] += er * _[i - 1];
};
auto polyval = [&](vector<double> const &c, double const &_x) -> double {
double res = 0.0;
int m = sz(c);
for (int ii = 0; ii < m; ii++)
res += c[ii] * pow(_x, (double)(m - ii - 1));
return res;
};
int __ = 1;
//input(_);
while (__--)
{
double _x, temp;
cin >> _x;
vector<double> x, y;
while (cin >> temp)
x.emplace_back(temp), y.emplace_back(cos(temp));
n = x.size();
vector<vector<double>> a(n, vector<double>(n));
int i, j;
for0(i, n) a[i][0] = y[i];
forab(j, 1, n - 1) forab(i, j, n - 1) a[i][j] = (a[i][j - 1] - a[i - 1][j - 1]) / (x[i] - x[i - j]);
vector<double> v;
v.emplace_back(a[n - 1][n - 1]);
forba(i, 0, n - 2)
{
polymul(v, -x[i]);
int l = sz(v);
v[l - 1] += a[i][i];
}
for0(i, n)
pr("%g\n", v[i]);
for0(i, n)
{
for0(j, n)
pr("%g ", a[i][j]);
puts("");
}
double _y = polyval(v, _x);
pr("%g\n", _y);
pr("%.6e",fabs(_y-cos(_x)));
}
return 0;
}
python代码
'''
Author: csc
Date: 2021-04-29 23:00:57
LastEditTime: 2021-04-30 09:58:07
LastEditors: Please set LastEditors
Description: In User Settings Edit
FilePath: \code_py\newton_.py
'''
import numpy as np
def difference_quotient(x, y):
n = len(x)
a = np.zeros([n, n], dtype=float)
for i in range(n):
a[i][0] = y[i]
for j in range(1, n):
for i in range(j, n):
a[i][j] = (a[i][j-1]-a[i-1][j-1])/(x[i]-x[i-j])
return a
def newton(x, y, _x):
a = difference_quotient(x, y)
n = len(x)
s = a[n-1][n-1]
j = n-2
while j >= 0:
s = np.polyadd(np.polymul(s, np.poly1d(
[x[j]], True)), np.poly1d([a[j][j]]))
j -= 1
for i in range(n):
print('%g' % s[n-1-i])
for i in range(n):
for j in range(n):
print('%g' % a[i][j], end=' ')
print()
_y = np.polyval(s, _x)
print('%g' % _y)
# re_err
real_y = np.cos(_x)
err = abs(_y-real_y)
print('%.6e' % err)
def main():
_x = float(input())
x = list(map(float, input().split()))
y = np.cos(x)
newton(x, y, _x)
if __name__ == '__main__':
main()
来源:https://blog.csdn.net/weixin_45720246/article/details/116300384
猜你喜欢
- ckptfrom tensorflow.python import pywrap_tensorflow checkpoint_path =
- 不止一次在微信、知乎有读者朋友跑过来问:看完了基础书,甚至看两遍了,但自己写的时候还是没思路,我该怎么办?编程在我看来就是一门手艺活,绝不是
- 均匀性度量图像分割是图像像素分割的一种方法,当然还有其他很多的方法。这里简单的介绍下其原理和实现代码【有源码】其流程大概分为一下几步1、确定
- poplib模块接收邮件python的poplib模块是用来从pop3收取邮件的,也可以说它是处理邮件的第一步。POP3协议并不复杂,它也是
- 看到了很多关于如何读出图片的长度的高度的方法,其实都不实用,大多都是通过图片的大小来判断的,图片的种类众多,通过大小来判断难免要制造很多的代
- 注意,在改变数值之前锁定应用,确保一段时间里只有一个客户执行该语句。<SCRIPT LANGUAGE="VBScr
- 问:如何让内联框架里面的网页背景透明?(用iframe嵌套一网页,怎么能够设置其背景为透明以便能显示父框架中网页的背景图?)答:这是需要 I
- 笔者日积月累了许多精彩、实用的Web特效的制作,这些特效几乎都是比较常用的网页特效。现在我就把这些经过
- 正在看的ORACLE教程是:Oracle数据库集复制方法浅议。前言 日益增长的分布式应用需求要求实现更好分布式的软件环境,不断推动着分布式技
- 不错,这个是一个文章详细页,没有左右两栏布局,不过这里我重点要讲的是合理的布局,在稍后的文章中我会详细的介绍浮动元素。好,回到刚才的话题,大
- 首先,我们先来看看,如果是人正常的行为,是如何获取网页内容的。(1)打开浏览器,输入URL,打开源网页(2)选取我们想要的内容,包括标题,作
- 将转储设备加入到SQL Server备份数据库的地方。在SEM中转储设备是可见性的,并且在设备上的信息被存储在主要数据库的sysdevice
- 列表生成式语法:[x*x for x in range(0,10)] //列表生成式,这里是中括号//结果 [0, 1, 4, 9, 16,
- 1、PHP中的抽象类PHP 5 支持抽象类和抽象方法。定义为抽象的类不能被实例化。任何一个类,如果它里面至少有一个方法是被声明为抽象的,那么
- 最近在研究WEB布局,遇到<H1>标签了,<H1>标签很重要。在一般教程中都这么说,<H1>标签在同一页
- 在Python中将字符串转换为集合使用 set() 类将字符串转换为集合,例如 my_set = set(my_str)。 set() 类将
- 本文将详细解释这些函数的使用方法。首先,我们介绍Python语言中类似于Windows系统的dir命令的列出文件功能,然后描述如何测试一个文
- 这是个郁闷的问题。主级获得ID列表 select ID from FS_SD_Address where PID=0
- 一、概述公司新购了一批PC,准备把几个性能较优的PC升级为数据库服务器,替换老旧的机器。公司有套POS终端软件,后台数据存储是 MySQL
- URL 编码是什么东东呢?看看我从网上抄的定义: 引用: url编码是一种浏览器用来打包