Несколько проблем в вашем коде:
Во-первых, да, у вас нет отступа (но я предполагаю, что это из-за того, что оно не копируется должным образом, так как это приведет к ошибке, а не к неправильному значению).В будущем убедитесь, что отступ в вашем вопросе соответствует тому, что у вас есть на вашем компьютере, прежде чем публиковать ...
Тогда термин должен быть добавлен в for тогда и только тогда, когда он находится в соответствующемсумма ... Здесь вы помещаете все в двойной цикл for, что соответствует наличию всех членов в двойной сумме.
Наконец, range(1,n) уже останавливается только на n-1, поэтому вы хотите удалить те-1 в диапазонах.
В конце концов:
def double_integral(f,a,b,c,d,nx,ny):
hx = (b-a)/nx
hy = (d-c)/ny
first_term = (f(a,c)+f(a,d)+f(b,c)+f(b,d))
i_sum = 0
for i in range(1,ny):
i_sum += f(a,c+i*hy)+f(b, c+i*hy)
j_sum = 0
for j in range(1,nx):
j_sum += f(a+j*hx,c)+f(a+j*hx,d)
ij_sum = 0
for i in range(1,ny):
for j in range(1,nx):
ij_sum += f(a+j*hx,c+i*hy)
integral = (first_term/4 + i_sum/2 + j_sum/2 + ij_sum) * hx * hy
return integral
def t(x,y):
return x*(y**(2))
print(double_integral(t,0,2,0,1,10,10))
0.6700000000000003
Вы приблизитесь к 2/3, выбрав количество шагов больше, чем 10 ...
И вы можете быть более питоническим, используя понимание суммы:
def double_integral(f,a,b,c,d,nx,ny):
hx = (b-a)/nx
hy = (d-c)/ny
first_term = (f(a,c)+f(a,d)+f(b,c)+f(b,d))
i_sum = sum(f(a,c+i*hy)+f(b, c+i*hy) for i in range (1,ny))
j_sum = sum(f(a+j*hx,c)+f(a+j*hx,d) for j in range(1,nx))
ij_sum = sum(f(a+j*hx,c+i*hy) for i in range (1,ny) for j in range(1,nx))
integral = (first_term/4 + i_sum/2 + j_sum/2 + ij_sum) * hx * hy
return integral