Проблема с числовой интеграцией

В настоящее время я застрял на действительно простой ошибке, спрятанной где-то, и, надеюсь, кто-то сможет пролить свет на это.

Я пытаюсь численно интегрировать набор данных с обычными методами, и получаю результаты, которые я не совсем ожидал.

Затем я вернулся к первым принципам и создал следующий код. По некоторым причинам я не могу произвести 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] станет бесконечно малым, согласно теоретическому определению интеграла

Основываясь на этом определении интеграла,

\ int_a ^ b! Р (х) ах

где F - антипроизводный или неопределенный интеграл от f.

Для тебя y3 или же y4, в x=(x[i] + x[i+1])/2 у вас есть у:

f ((x [i] + x [i + 1]) / 2) ~ = F (x [i + 1]) - F (x [i])

Для тебя y2, в x=(x[i] + x[i+1])/2 у вас есть у:

f ((x [i] + x [i + 1]) / 2) = F ((x [i] + x [i + 1]) / 2)

Когда "x[i+1] - x[i] становится бесконечно малым", вы должны иметь

F ((x [i] + x [i + 1]) / 2) ~ = F (x [i]) ~ = F (x [i + 1])

F (x [i + 1]) - F (x [i]) ~ = 0

Другими словами, когда "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,

Другие вопросы по тегам