詳解Python牛頓插值法
一、牛頓多項(xiàng)式
拉格朗日多項(xiàng)式的公式不具備遞推性,每個(gè)多項(xiàng)式需要單獨(dú)構(gòu)造。但很多時(shí)候我們需要從若干個(gè)逼近多項(xiàng)式選擇一個(gè)。這個(gè)時(shí)候我們就需要一個(gè)具有遞推關(guān)系的方法來(lái)構(gòu)造——牛頓多項(xiàng)式
這里的的a0,a1…等可以通過(guò)逐一帶入點(diǎn)的值來(lái)求得。但是當(dāng)項(xiàng)數(shù)多起來(lái)時(shí),會(huì)發(fā)現(xiàn)式子變得很大,這個(gè)時(shí)候我們便要引入差商的概念(利用差分思想)具體見(jiàn)式子能更好理解
這里在編程實(shí)現(xiàn)中我們可以推出相應(yīng)的差商推導(dǎo)方程
d(k,0)=y(k)
d(k,j)=(d(k,j-1)-d(k-1,j-1)) / (x(k)-x(k-j))
二、例題
【問(wèn)題描述】考慮[0,3]內(nèi)的函數(shù)y=f(x)=cos(x)。利用多個(gè)(最多為6個(gè))節(jié)點(diǎn)構(gòu)造牛頓插值多項(xiàng)式。
【輸入形式】在屏幕上依次輸入在區(qū)間[0,3]內(nèi)的一個(gè)值x*,構(gòu)造插值多項(xiàng)式后求其P(x*)值,和多個(gè)節(jié)點(diǎn)的x坐標(biāo)。
【輸出形式】輸出牛頓插值多項(xiàng)式系數(shù)向量,差商矩陣,P(x*)值(保留6位有效數(shù)字),和與真實(shí)值的絕對(duì)誤差(使用科學(xué)計(jì)數(shù)法,保留小數(shù)點(diǎn)后6位有數(shù)字)。
【樣例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說(shuō)明】
輸入:x為0.8,3個(gè)節(jié)點(diǎn)為(k, cos(k)),其中k = 0, 0.5, 1。
輸出:
牛頓插值多項(xiàng)式系數(shù)向量,表示P2(x)=-0.429726x^2 - 0.0299721x + 1;
3行3列的差商矩陣;
當(dāng)x為0.8時(shí),P2(0.8)值為0.700998
與真實(shí)值的絕對(duì)誤差為:4.291237*10^(-3)
【評(píng)分標(biāo)準(zhǔn)】根據(jù)輸入得到的輸出準(zhǔn)確
三、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()
到此這篇關(guān)于詳解Python牛頓插值法的文章就介紹到這了,更多相關(guān)Python牛頓插值法內(nèi)容請(qǐng)搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家!
相關(guān)文章
python用for循環(huán)求和的方法總結(jié)
在本篇文章里小編給各位分享了關(guān)于python用for循環(huán)求和的方法以及相關(guān)實(shí)例代碼,需要的朋友們參考學(xué)習(xí)下。2019-07-07python實(shí)現(xiàn)顏色rgb和hex相互轉(zhuǎn)換的函數(shù)
這篇文章主要介紹了python實(shí)現(xiàn)顏色rgb和hex相互轉(zhuǎn)換的函數(shù),可實(shí)現(xiàn)將rgb表示的顏色轉(zhuǎn)換成hex值的功能,非常具有實(shí)用價(jià)值,需要的朋友可以參考下2015-03-03python實(shí)現(xiàn)簡(jiǎn)易云音樂(lè)播放器
這篇文章主要介紹了python實(shí)現(xiàn)簡(jiǎn)易云音樂(lè)播放器,具有一定的參考價(jià)值,感興趣的小伙伴們可以參考一下2018-01-01Django開(kāi)發(fā)web后端對(duì)比SpringBoot示例分析
這篇文章主要介紹了Django開(kāi)發(fā)web后端對(duì)比SpringBoot示例分析,有需要的朋友可以借鑒參考下,希望能夠有所幫助,祝大家多多進(jìn)步,早日升職加薪2023-12-12matplotlib實(shí)現(xiàn)矩陣和圖像的可視化表示
這篇文章主要為大家詳細(xì)介紹了如何利用matplotlib實(shí)現(xiàn)矩陣和圖像的可視化表示,文中的示例代碼講解詳細(xì),具有一定的學(xué)習(xí)價(jià)值,感興趣的小伙伴可以了解下2024-03-03Tensorflow2.1 完成權(quán)重或模型的保存和加載
這篇文章主要為大家介紹了Tensorflow2.1 完成權(quán)重或模型的保存和加載,有需要的朋友可以借鑒參考下,希望能夠有所幫助,祝大家多多進(jìn)步,早日升職加薪2022-11-11