概念
设f(x)在N+1个点(x0,y0)…(xn,yn)处的值已知,其中值xk在区间【a,b】上,xk互不相同,满足a≤x0<x1<…<xn≤b,yk=f(xk)。求任一插值点x对应的y。
插值法的思路在于通过构造一个N次多项式P(x),使其通过N+1个点,然后求x处的函数值y。
拉格朗日逼近
例题
【问题描述】考虑[0.0,1.2]内的函数y=f(x)=cos(x)。利用多个(2,3,4等)节点构造拉格朗日插值多项式。
【输入形式】在屏幕上依次输入在区间[0.0,1.2]内的一个值x*,构造插值多项式后求其P(x*)值,和多个节点的x坐标。
【输出形式】输出插值多项式系数矩阵,拉格朗日系数多项式矩阵和P(x*)值(保留6位有效数字)。
【样例1输入】
0.3
0 0.6 1.2
【样例1输出】
-0.400435
-0.0508461
1
1.38889 -2.5 1
-2.77778 3.33333 -0
1.38889 -0.833333 0
0.948707
【样例1说明】输入:x为0.3,3个节点的x坐标分别为x0=0,x1=0.6和x2=1.2。
输出:插值多项式系数矩阵,则插值多项式P2(x)表示为-0.400435x**2-0.0508461x+1;
拉格朗日系数多项式矩阵,则P2(x)表示为:y0(1.38889x2-2.5x+1)+y1*(-2.77778x2+3.33333x-0)+y2*(1.38889x**2-0.83333x+0);
当x*为0.3时,P2(0.3)值为0.948707。
【评分标准】根据输入得到的输出准确
ACcode:
c++ code(后面有python代码)
#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 pdd pair<double, double>#define vd vector<double>#define vdd vector<vector<double>>#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;vd operator*(vd const &y, vdd const &l){int m = sz(y);vd c(m);for (int ii = 0; ii < m; ii++)for (int jj = 0; jj < m; jj++)c[ii] += y[jj] * l[jj][ii];return c;}signed main(){int __ = 1;//input(__);auto polymul = [&](vd &v, double const &d, double const &err) {int m = sz(v);vd _ = v;_.emplace_back(0);for (int g = 0; g < m; g++)_[g + 1] += v[g] * d, _[g] /= err;_[m] /= err;v = _;};auto polyval = [&](vd 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;};while (__--){double _x;cin >> _x;vd x, y;double X, Y;while (cin >> X)Y = cos(X), x.emplace_back(X), y.emplace_back(Y);int n = sz(x);int i, j, k;vdd l(n, vd(n));for0(k, n){vd v(1, 1);for0(j, n){if (j != k){double d = -x[j];double err = x[k] - x[j];polymul(v, d, err);}}l[k] = v;}vd c = y*l;for0(i, n)pr("%g\n", c[i]);for0(i, n){for0(j, n)pr("%g ", l[i][j]);puts("");}pr("%g\n", polyval(c, _x));}return 0;}
python --code
'''Author: cscDate: -04-24 10:52:25LastEditTime: -04-24 14:06:57LastEditors: Please set LastEditorsDescription: In User Settings EditFilePath: \code_py\lagrange.py'''import numpy as npdef lagrange(_x, n, x, y):l = np.zeros([n, n], dtype=np.double)for k in range(n):v = 1for j in range(n):if j != k:v = np.polymul(v, np.poly1d([1, -x[j]]))/(x[k]-x[j])l[k, :] = vc = y@lfor i in range(n):print('%g' % c[i])for i in range(n):for j in range(n):print('%g' % l[i][j], end=' ')print()_y = np.polyval(c, _x)print('%g' % _y, end=' ')def main():_x = float(input())x = np.array(list(map(float, input().split())))n = len(x)y = np.cos(x)lagrange(_x, n, x, y)if __name__ == '__main__':main()