In [2]:
import numpy as np
from matplotlib import pylab as plt
from scipy.integrate import trapz
%matplotlib inline
In [3]:
def rect_int(f, x):
h = x[1] - x[0]
a = x[0]
b = x[-1]
x_mids = np.linspace(a+h/2, b-h/2, x.size-1)
return h*np.sum(f(x_mids))
def f(x):
return x**2
In [4]:
a = 0
b = 5
n = 5
x = np.linspace(a,b,n)
y = f(x)
plt.plot(x,y)
print(rect_int(f, x))
print(np.trapz(f(x),x))
In [5]:
x[:-1]
Out[5]:
In [6]:
x
Out[6]:
In [7]:
h = 1.25
np.linspace(a+h/2, b-h/2, x.size-1)
Out[7]:
In [8]:
def g(x):
x1 = x[0]
x2 = x[1]
return x1 + x2**2
def grid_vect(a,b,N):
x = np.linspace(a, b, N, endpoint=True)
h = x[1] - x[0]
return np.mgrid[a:b+h:h, a:b+h:h].reshape(2,-1).T
def double_rect_int(f, x):
temp = x[0]
#a = temp[0]
#b = temp[-1]
h = temp[1] - temp[0]
# x = x[:,:-N]
# x1 = x[0]
# x2 = x[1]
# x_grid = np.meshgrid(x1,x2)
# z = f(x)
return res
xy = np.transpose(grid_vect(0,5,5))
In [9]:
g(xy)
Out[9]:
In [10]:
xy[:,:-5]
Out[10]:
In [11]:
xy
plt.scatter(xy[0],xy[1])
Out[11]:
In [12]:
x_grid = np.meshgrid(x,x)
In [13]:
x2g = x_grid[1]
x1g = x_grid[0]
In [14]:
x1g
Out[14]:
In [15]:
x_grid
Out[15]:
In [16]:
x1g.shape
Out[16]:
In [22]:
print(x_grid)
print('========')
fff1 = x_grid[0].reshape(n**2)
fff2 = x_grid[0].reshape(n**2)
print(fff1)
print(fff2)
print('========')
print(fff1.reshape((5,5)))
In [18]:
plt.scatter(x1g,x2g)
Out[18]:
In [ ]:
In [19]:
from mpl_toolkits.mplot3d.axes3d import Axes3D
fig = plt.figure(figsize=(14,6))
ax = fig.add_subplot(1, 2, 1, projection='3d')
p = ax.plot_surface(x_grid[0], x_grid[1], g(x_grid), rstride=4, cstride=4, linewidth=0)
In [20]:
xy
Out[20]:
In [21]:
xy_y = xy[1]
print(xy_y)
In [24]:
x_grid
Out[24]:
In [23]:
g(x_grid)
Out[23]:
In [28]:
xy.reshape((2,5,5))
Out[28]:
In [26]:
g(xy) - (g(x_grid)).reshape(n**2)
Out[26]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [97]:
def rectangular(f, a, b, n):
h = float(b - a)/n
x = np.linspace(a+h/2, b-h/2, n)
return h*np.sum(f(x))
def rectangular_double2(f, a, b, c, d, nx, ny):
"""
Вычисляет значение двойного интеграла по формуле прямоугольников
с использованием функции для одномерного интеграла
f - подынтегральная функция
a, b, c, d - границы прямоугольной области интегрирования
nx, ny - количество частичных отрезков по x и y соответственно
"""
g = lambda x: rectangular(lambda y: f(x, y), c, d, ny)
return rectangular(g, a, b, nx)
aaa = lambda x,y: x**2+y**3
rectangular_double2(aaa, -1,1,-1,1,50,50)
Out[97]:
In [96]:
#bbb =
NN = 50
xt = np.linspace(-1,1,NN)
yt = np.linspace(-1,1,NN)
xyt = np.transpose(grid_vect(-1,1,NN))
Zt = aaa(xyt[0], xyt[1])
Zt = Zt.reshape((NN,NN))
#print(xyt.shape, Zt.shape)
a1 = np.trapz(Zt,xyt[1].reshape((NN,NN)))
#a1 = np.trapz(Zt, xyt, axis = -1)
#print(a1)
a2 = np.trapz(a1, xt)
print(a2)
#print(Zt)
#print(xyt[0])
#print(aaa(xyt[0],xyt[1]).reshape(NN**2))
#print(xyt)
In [65]:
xyt
Out[65]:
In [ ]:
In [ ]:
In [ ]:
In [47]:
np.transpose(np.array([xt,yt]))
Out[47]:
In [ ]:
In [ ]:
In [ ]:
In [32]:
def rectangular(f, a, b, n):
h = float(b - a)/n
x = np.linspace(a+h/2, b-h/2, n)
return h*np.sum(f(x))
def rectangular_double2(f, a, b, c, d, nx, ny):
g = lambda x: rectangular(lambda y: f(np.array([[x],[y]])), c, d, ny)
return rectangular(g, a, b, nx)
aaa = lambda XX: XX[0]**2+XX[1]**2
rectangular_double2(aaa, -1,1,-1,1,50,50)
Out[32]: