Ссылка на проблему: 1
\ int_b ^ a f_1 (x) e ^ {\ int_c ^ x f_2 (y) dy} dx.
Это легко сделать с помощью Mathematica, а также существуют решения в Matlab. Но я не могу понять, как оценить этот интеграл с помощью scipy.quad или аналогичных функций. Поэтому я попытался реализовать это с помощью составного правила Симпсона, но в некоторых случаях ошибка велика. Код адаптирован из численного анализа Бремена и Фэйрса. Например, f1 (x) = co sh (x) и f2 (x) = 2x + sinh (x), a = 0, b = 4, c = 4. Из ММА я получаю 0,7699, в то время как с этой числовой схемой я получаю 0,759206.
def compSimpsInt(fun1, fun2, a, b, c):
"""
Numerical integration of the form
\int_a^b fun1(x,y) e^{\int_c^x fun2(y)dy} dx
"""
n = len(fun1)
m = len(fun2)
h = (b - a) / n
J1 = 0
J2 = 0
J3 = 0
for i in range(n):
x = a + (i) * h
HX = (x - c) / m
K1 = fun2[-1] + fun2[i] # end terms. fun[-1]: at L, fun2[i] at x.
K2 = 0 # Even terms
K3 = 0 # Odd terms
Q = np.interp(np.linspace(c,x,m), np.linspace(a,b,n), fun2)
K2 = sum(Q[i] for i in range(1, m-1) if i % 2 != 0)
K3 = sum(Q[i] for i in range(1, m-1) if i % 2 == 0)
L = (K1 + 2 * K2 + 4 * K3) * HX / 3
LL = np.exp(L) * fun1[i]
if i == 0 or i == n:
J1 += LL
elif i % 2 == 0:
J2 += LL
else:
J3 += LL
J = (J1 + 2 * J2 + 4 * J3) * h / 3
return J
Буду признателен, если кто-нибудь поможет, либо укажет на ошибку в коде, либо использует scipy.integrate (или любой другой пакет). ), если интеграция может быть аппроксимирована ..