Разложить функцию в ряд

python

numpy, scipy — нужные тебе кейворды для гугла и stackoverflow

Общий ответ дать не могу, но вот родственная задача:

Даны функции f, f1 … fm и сетка x[1] … x[n]. Найти оптимальные коэффициенты a1 … am такие, чтобы остаток a1 * f(x[i]) + ... am * fm(x[i]) - f(x[i]) был близок к нулю для каждого x[i] на сетке.

Задачу можно решить в смысле наименьших квадратов, если посчитать каждую функцию в каждой точке, посчитать искомую f(x[i]) в каждой точке сетки, после чего передать получившуюся матрицу и получившийся вектор в numpy.linalg.lstsq.

Вашу задачу можно свести к ней, если задать сетку и слегка изменить функции (смещением аргумента на константу).

Задачу можно решить в смысле наименьших квадратов, если посчитать каждую функцию в каждой точке, посчитать искомую f(x[i]) в каждой точке сетки, после чего передать получившуюся матрицу и получившийся вектор в numpy.linalg.lstsq.

Почему есть решение? Почему единственное?

Гуглите многокомпонентный ряд Тейлора, Вам нужно разложение до линейного члена. Там даже чуть проще будет: f0 + a1* x1+…+an* fn.

Что бы его построить надо считать производные, это умеет например sympy

Почему есть решение? Почему единственное?

Ерунду спрашиваю, документация:

Computes the vector x that approximately solves the equation a @ x = b. The equation may be under-, well-, or over-determined (i.e., the number of linearly independent rows of a can be less than, equal to, or greater than its number of linearly independent columns). If a is square and of full rank, then x (but for round-off error) is the “exact” solution of the equation. Else, x minimizes the Euclidean 2-norm ||b - ax||. If there are multiple minimizing solutions, the one with the smallest 2-norm ||x|| is returned.

ЧЯДНТ?

import numpy
from scipy import linalg
from matplotlib import pyplot

def foo(t):
    return t
def res(t, x, n, f):
    res = 0.0
    for i in range(n):
        for j in range(len(f)):
            res += x[i * len(f) + j] * f[j](t)
    return res

def f0(t):
    return 1
def f1(t):
    return t
def f2(t):
    return t ** 2

a = 0
b = 5
f = [f0, f1, f2]
n = 5
m = 5
n_columns = n * len(f)
mat = numpy.eye(m, n_columns, dtype=float)
vec = numpy.arange(m, dtype=float)
for i in range(m):
    xi = i * (b - a) / m
    vec[i] = foo(xi)
    for j in range(n):
        dx = j * (b - a) / n
        for k in range(len(f)):
            mat[i, j * len(f) + k] = f[k](xi + dx)
# mat * x = vec
x, residuals, rank, s = linalg.lstsq(mat, vec)

t = numpy.linspace(a, b, m)
pyplot.plot(t, foo(t), color='red')
pyplot.plot(t, res(t, x, n, f), color='blue')
pyplot.show()

ЧЯДНТ?

В Вашем случае получилась неполноранговая матрица, на ней эти методы плохо работают.

n_columns = n * len(f)

А почему получилось несколько столбцов на функцию? Я думал, на каждую функцию f[i] получается по одному смещению x[i]:

import numpy as np
import matplotlib.pyplot as plt

# foo(x) ~ a1 * f1(x + x1) + … + an * fn(x + xn)

def foo(x): return x

# f[i](x + dx[i]) как выше
f = [
	lambda x: np.ones_like(x),
	lambda x: x,
	lambda x: x ** 2,
]
dx = [.3, .2, .1]
assert(len(f) == len(dx))

# целевой диапазон
a, b = 0, 5

# вручную заданная сетка заданного размера
m = 5
x = np.linspace(a, b, m)

# вычисляем матрицу
A = np.empty((m, len(f)))
for i in range(len(f)): A[:,i] = f[i](x + dx[i])

coef, resid, _, _ = np.linalg.lstsq(A, foo(x), rcond = None)

print(coef)

plt.plot(x, foo(x), label = 'orig')
plt.plot(x, A @ coef, label = 'pred')
plt.legend()
plt.title('2-norm of residuals = %g' % resid)
plt.show()

Нашел ошибку, не учитывал сдвиг:

def res(t, x, n, f):
    res = 0.0
    for i in range(n):
        dx = i * (b - a) / n
        for j in range(len(f)):
            res += x[i * len(f) + j] * f[j](t + dx)
    return res

А почему получилось несколько столбцов на функцию?

Из каких соображений выбирать dx? Я раскладываю по всем f[i](x + j * (b - a) / n)

Я раскладываю по всем f[i](x + j * (b - a) / n)

Что, разумеется, неверно в изначальной формулировке.

Если разложение не допускает совпадений fi(x + dx1) + fi(x + dx2) то
для каждого набора dx решать систему и искать минимальную невязку?
Или можно как-то упростить задачу?

Из каких соображений выбирать dx?

А для чего всё это делается? Может, и линейная система тут ни к чему, и задачу нужно решать совсем по-другому? Здесь-то, понятно, самое простое решение будет с dx = 0, но в общем случае это не так. Никогда не занимался этим вопросом, но, может, нужно заняться символьной регрессией? Найти такие dx, при которых foo(x) линейна по каждой из функций f[i](x + dx[i]). Но что делать, если их нет?

Если разложение не допускает совпадений fi(x + dx1) + fi(x + dx2) то для каждого набора dx решать систему и искать минимальную невязку?

Я ошибся, когда сказал, что дело в недостаточном ранге матрицы предикторов. Это, в принципе, плохо (делает решение неоднозначным), но numpy.linalg.lstsq с этим справляется, выбирая решение с минимальной нормой. Так что даже если fi(x + dx1) и fi(x + dx2) - это не только одинаковые функции с разным dx, но и линейно зависимые друг от друга (fi(x + dx1) * const == fi(x + dx2)), какие-то коэффициенты Вам всё равно достанутся, и невязку они будут минимизировать. Серьёзной проблемы здесь нет.

для каждого набора dx решать систему и искать минимальную невязку?

Формально, можно считать это задачей оптимизации. Записать невязку ||A(f,x,dx) @ coef - foo(x)|| как функцию от вектора dx (остальные параметры зафиксированы заранее или вычисляются из них) и в таком виде запихнуть в метод минимизации, авось что-нибудь скажет. Нужно готовиться к тому, что ответов с примерно одинаковой невязкой будет много (например, в примере про foo(t) = t через 1, t и t^2, пока ошибки округления от работы с большим t^2 всё не сломают, разницы от выбора dx практически нет).

А для чего всё это делается?

Для развлечения (:

Спасибо за ответы!

Для этого никакие программы не нужны: для любых f, f_1…f_n, x_1…x_n

foo(f, f_1…f_n)=(0,x_1)….(0,x_n)

f(x)=eps(x)

Вот и готовый ответ.