Полином, удовлетворяющий интегралу и двум точкам - PullRequest
Новая
       86

Полином, удовлетворяющий интегралу и двум точкам

1 голос
/ 24 марта 2020 
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, но у меня также есть проблемы с реализацией.

Ответы [ 2 ]

3 голосов
/ 24 марта 2020

Используя математическую библиотеку sympy , Python symboli c, проблему можно сформулировать следующим образом: 

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 уравнениями. Уравнения должны быть «вычислены вручную», что более подвержено ошибкам и более тщательно разработано, чтобы распространяться на более высокие степени. 

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)
1 голос
/ 24 марта 2020

Вы можете решить это в общем, воспользовавшись тем, что

$$ \int_{-1}^1 (a x^2 + b x + c) dx = 2/3 (a + 3 c) $$

и добавив его в качестве ограничения. Тогда у вас есть 3 уравнения для 3 неизвестных (a, b, c).

Есть и другие интересные приемы, это интересный вопрос. Попробуйте поиграть с написанием формулы в терминах a(x-b)(x-c), тогда у вас будет 3bc + 1 = 0. Также любое решение, начинающееся с точек (x0,y0),(x1,x1), имеет аналогичное решение для (k*x0,k*y0),(k*x1,k*y1).

...