2000字范文,分享全网优秀范文,学习好帮手!
2000字范文 > 【计算方法】拉格朗日插值法

【计算方法】拉格朗日插值法

时间:2019-11-02 17:09:18

相关推荐

【计算方法】拉格朗日插值法

概念

设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()

本内容不代表本网观点和政治立场,如有侵犯你的权益请联系我们处理。
网友评论
网友评论仅供其表达个人看法,并不表明网站立场。