Проблемы с простым численным интегрированием в Python ⇐ Python

[Anonymous]

Я пытаюсь решить простой интеграл на Python, используя различные методы численного интегрирования. Однако когда я пытаюсь вычислить интегралы, я получаю неправильные результаты для малых значений моей координаты x. Это путаница, потому что я пытаюсь решить некоторую систему уравнений, в которой одна из переменных, которые я эволюционирую во времени, является интегралом некоторых величин. Конечно, у меня нет функциональной формы для выражения, которое я хочу интегрировать, но вместо этого у меня есть массив данных (по одному f для каждого x). Знаете ли вы, как получить более точные результаты моих интегралов?
Здесь я представляю простой пример проблемы, с которой я столкнулся. Предположим, у меня есть некоторые данные, которые на данный момент я буду рассматривать только простой x^7 (вообще у меня нет функциональной формы для функций, которые я хочу интегрировать, и это всего лишь пример), и предположим, что я хочу интегрировать такую ​​функцию. Аналитическое решение интеграла простое: Int = x^8/8. Затем, если я напечатаю или построю график 8*Int/x^8, у меня должно быть то, что каждая запись вывода должна быть равна единице. Затем позвольте мне рассмотреть следующий код:

Код: Выделить всё

import numpy as np
from scipy.integrate import simps

x = np.linspace(0,10,101)
f = x**7
Int = np.array([simps(f[:i+1], x[:i+1]) for i in range(len(x))])
Div = 8*Int/(x**8)
Div[0] = 1 #To avoid division over zero
print(Div)

Вывод этого простого кода:

Код: Выделить всё

[1.           4.           1.375        1.2687090383 1.033203125
1.05997312   1.0069158665 1.0195120352 1.0022277832 1.0080582909
1.00092      1.0038955372 1.0004456376 1.0021034429 1.0002411835
1.0012324471 1.0001416206 1.0007688345 1.0000885172 1.0005040075
1.000058125  1.0003439529 1.0000397249 1.0002426572 1.0000280617
1.0001760417 1.000020381  1.0001307892 1.000015157  1.0000991831
1.0000115043 1.000076571  1.0000088885 1.0000600497 1.0000069756
1.0000477527 1.0000055507 1.0000384478 1.0000044717 1.0000313025
1.0000036426 1.0000257424 1.000002997  1.0000213637 1.0000024883
1.0000178775 1.0000020831 1.0000150742 1.0000017571 1.0000127994
1.0000014925 1.0000109379 1.0000012758 1.0000094029 1.0000010971
1.0000081279 1.0000009486 1.0000070618 1.0000008244 1.000006165
1.0000007199 1.000005406  1.0000006314 1.0000047604 1.0000005561
1.0000042083 1.0000004917 1.000003734  1.0000004364 1.0000033246
1.0000003886 1.0000029699 1.0000003472 1.0000026612 1.0000003112
1.0000023916 1.0000002797 1.0000021553 1.0000002521 1.0000019475
1.0000002278 1.0000017641 1.0000002064 1.0000016018 1.0000001874
1.0000014577 1.0000001706 1.0000013295 1.0000001556 1.0000012151
1.0000001422 1.0000011127 1.0000001303 1.0000010209 1.0000001195
1.0000009384 1.0000001099 1.000000864  1.0000001012 1.0000007969
1.0000000933]

Как видите, большинство результатов очень близки к 1 (как и ожидалось), за исключением точек, которые очень близки к точке x = 0. Мне нужно иметь хорошее контроль над этими интегралами, поскольку я использую их для численного решения системы уравнений. Проблема с такими плохими результатами заключается в том, что в моей числовой схеме они дают мне некоторые градиенты, которые разрушают мои численные решения. Итак, вы знаете, как с этим справиться?
PS: Если вместо использования simps я использую mpmath, я могу получить правильный результат. Однако mpmath работает только с уже известными нам функциями (как в этом простом примере функции x^7). Однако, как я уже упоминал ранее, на самом деле я не знаю функциональной формы выражения, которую не хочу интегрировать, и вместо этого у меня есть только массив данных.
PS2: Я также можно попытаться увеличить размер моего массива x. Но даже если я это сделаю, у меня всегда будут ошибки в первых точках, ближайших к x = 0, и эти ошибки в конечном итоге испортят мои численные решения.
С уважением, LP

Подробнее здесь: https://stackoverflow.com/questions/786 ... -in-python