Проблема с числовой интеграцией
В настоящее время я застрял на действительно простой ошибке, спрятанной где-то, и, надеюсь, кто-то сможет пролить свет на это.
Я пытаюсь численно интегрировать набор данных с обычными методами, и получаю результаты, которые я не совсем ожидал.
Затем я вернулся к первым принципам и создал следующий код. По некоторым причинам я не могу произвести y3 или y4, чтобы соответствовать y2.
import numpy as np
import matplotlib.pyplot as plt
import math
f1 = 'x**2 - 1' # Starting function
start, stop, step = -2, 2, 21
x = np.linspace(start, stop, step)
xdash = []
xDoubleDash = []
y1 = eval(f1)
y2 = [] # integrated function
y3 = [] # box rule
y4 = [] # trapezium rule
y5 = [] # numerical differentialtion of integrated function
for i in range(0,len(x)-1):
xdashElement = (x[i] + x[i+1])/2
y2.append((math.pow(xdashElement,3)/3) - xdashElement) # Integrated Function: (x**3)/3 - x
xdash.append(xdashElement)
y3.append((y1[i]+y1[i+1])/2 * abs(((x[i] - x[i+1])))) # box rule
y4.append(((y1[i]+y1[i+1]) * abs(x[i] - x[i+1]))/2) # trapezium rule
for i in range(0,len(y2)-1):
xDoubleDashElement = (xdash[i] + xdash[i+1])/2
xDoubleDash.append(xDoubleDashElement)
y5.append((y2[i+1]-y2[i])/(xdash[i+1] - xdash[i]))
plt.plot(x, y1, 'b-')
plt.plot(xdash, y2, 'r-')
plt.plot(xdash, y3, 'g-')
plt.plot(xdash, y4, 'm-')
plt.plot(xDoubleDash, y5, 'c-')
plt.grid()
plt.show()
Заранее спасибо.
РЕДАКТИРОВАТЬ: я добавил y5, численное дифференцирование y2 в качестве проверки вменяемости.
1 ответ
Ваш y2
является неопределенным интегралом, в то время как ваш y3
а также y4
являются определенным интегралом по серии интервалов. Я не ожидаю, что они соответствуют друг другу. Если вы хотите проверить свою числовую интеграцию с анти-производным, ваш y3
а также y4
должен соответствовать этому y2
:
ya = (math.pow(x[i],3)/3) - x[i]
yb = (math.pow(x[i+1],3)/3) - x[i+1]
y2.append(yb-ya) # Integrated Function: (x**3)/3 - x
Редактировать в ответ на комментарии по интегралу:
Я ожидал, что значения y3 и y4 будут очень похожи на значения y2 и будут теоретически совпадать, если приращение x[i+1] - x[i] станет бесконечно малым, согласно теоретическому определению интеграла
Основываясь на этом определении интеграла,
где F - антипроизводный или неопределенный интеграл от f.
Для тебя y3
или же y4
, в x=(x[i] + x[i+1])/2
у вас есть у:
Для тебя y2
, в x=(x[i] + x[i+1])/2
у вас есть у:
Когда "x[i+1] - x[i] становится бесконечно малым", вы должны иметь
Другими словами, когда "x[i+1] - x[i] становится бесконечно малым", ваш y3
или же y4
почти горизонтальная линия в y=0
, Почему ваш y2
соответствует линии в y=0
?
Если вы действительно хотите увидеть свой y3
а также y4
соответствовать форме y2
вам понадобится y0
или " произвольная постоянная интеграции", чтобы начать и сделать y3
а также y4
как это:
y0 = -.7
y3 = [y0] # box rule
y4 = [y0] # trapezium rule
for i in range(0,len(x)-1):
xdashElement = (x[i] + x[i+1])/2
y3.append((y1[i]+y1[i+1])/2 * abs(((x[i] - x[i+1]))) + y3[len(y3) - 1]) # box rule
y4.append(((y1[i]+y1[i+1]) * abs(x[i] - x[i+1]))/2 + y4[len(y4) - 1]) # trapezium rule
y3.pop(0)
y4.pop(0)
Затем с увеличением вашего step
или "приращение x[i+1] - x[i] становится бесконечно малым", эти y3
а также y4
будет соответствовать y2
,