满足积分和两点的多项式

问题描述 投票:1回答:2
Consider two points (x_0, f_0) and (x_1, f_1) 
let p(x) be the degree two polynomial for which
    p(x_0) = f_0
    p(x_1) = f_1
and the integral of p(x) from -1 to 1 is equal to 0

Write a function which accepts two arguments
    1. a length 2 NumPy vector 'x' of floating point values, with 'x[i]' containing the value of x_i,
    2. a length 2 NumPy vector 'f' of floating point values, with 'f[i]' containing the value of f_i,
and which returns
    a length 3 NumPy vector of floating point values containing the power series coefficients, in order from the highest order term to the constant term, for p(x)

我不确定从哪里开始。我最初的想法是让微分方程P(1)= P(-1)的初始值为p(x_0)= f_0和p(x_1)= f_1,但我在实现时也遇到了问题。

python math integration polynomials polynomial-math
2个回答
3
投票

使用Python的符号数学库sympy,该问题可以表述为:

from sympy import symbols Eq, solve, integrate

def give_coeff(x, f):
    a, b, c, X = symbols('a, b, c, X')
    F = a * X * X + b * X + c  # we have a second order polynomial
    sol = solve([Eq(integrate(F, (X, -1, 1)), 0),  # the integral should be zero (2/3*a + 2*c)
                 Eq(F.subs(X, x[0]), f[0]),        # filling in x[0] should give f[0]
                 Eq(F.subs(X, x[1]), f[1])],       # filling in x[1] should give f[1]
                (a, b, c))   # solve for a, b and c
    return sol[a].evalf(), sol[b].evalf(), sol[c].evalf()

import numpy as np

coeff = give_coeff(np.array([1, 2]), np.array([3, 4]))
print(coeff)

该代码甚至可以扩展为任意程度的多项式:

from sympy import Eq, solve, symbols, integrate


def give_coeff(x, f):
    assert len(x) == len(f), "x and f need to have the same length"
    degree = len(x)
    X = symbols('X')
    a = [symbols(f'a_{i}') for i in range(degree + 1)]
    F = 0
    for ai in a[::-1]:
        F = F * X + ai
    sol = solve([Eq(integrate(F, (X, -1, 1)), 0)] +
                [Eq(F.subs(X, xi), fi) for xi, fi in zip(x, f)],
                (*a,))
    # print(sol)
    # print(F.subs(sol).expand())
    return [sol[ai].evalf() for ai in a[::-1]]

import numpy as np

coeff = give_coeff(np.array([1, 2]), np.array([3, 4]))
print(coeff)
print(give_coeff(np.array([1, 2, 3, 4, 5]), np.array([3, 4, 6, 9, 1])))

PS:仅使用numpy求解二次方程,可以使用np.linalg.solve来求解带有3个方程的3个未知数的线性系统。这些方程需要“手工计算”,这更容易出错,并且要扩展到更高的程度,因此更加复杂。

np.linalg.solve

0
投票

您可以利用以下事实来一般性地解决此问题

$$ \ int _(-1)^ 1(a x ^ 2 + b x + c)dx = 2/3(a + 3 c)$$

并添加为约束。然后,您有3个方程式,涉及3个未知数(a,b,c)。

还有其他有趣的技巧,这是一个整洁的问题。尝试用import numpy as np def np_give_coeff(x, f): # general equation: F = a*X**2 + b*X + c # 3 equations: # integral (F, (X, -1, 1)) == 0 or (2/3*a + 2*c) == 0 # a*x[0]**2 + b*x[0] + c == f[0] # a*x[1]**2 + b*x[1] + c == f[1] A = np.array([[2/3, 0, 2], [x[0]**2, x[0], 1], [x[1]**2, x[1], 1]]) B = np.array([0, f[0], f[1]]) return np.linalg.solve(A, B) coeff = np_give_coeff(np.array([1, 2]), np.array([3, 4])) print(coeff) 来写公式,然后得到a(x-b)(x-c)。而且,任何以点3bc + 1 = 0开头的解都具有与(x0,y0),(x1,x1)类似的解。

© www.soinside.com 2019 - 2024. All rights reserved.