«Оптимал квадратур формулалар орқали Фурье интегралларини тақрибий ҳисоблаш алгоритми» дастурий тизими

'''
F0, F1 and F2 are Fourier transforms of
rect, linear and quadratic functions
'''

#h=0.1
w1 = np.arange(a,b+0.1,0.1)
M_1 = len(w1) - 1
h_1 = (b - a) / M_1
x1 = np.arange(a,b+0.1,0.1)
N1 = length(x1) - 1

C_1 = np.zeros(len(w1), len(x1))

E0_1 = np.zeros(1, len(w1))
E1_1 = np.zeros(1, len(w1))
E2_1 = np.zeros(1, len(w1))
E11_1 = np.zeros(1, len(w1))
E21_1 = np.zeros(1, len(w1))

f0_1 = np.zeros(1, len(x1))

f1_1 = np.zeros(1, len(x1))
f2_1 = np.zeros(1, len(x1))
f11_1 = np.zeros(1, len(x1))
f21_1 = np.zeros(1, len(x1))
for n in range(N1+1):
f0_1[n] = rectfunc(x1(n), a1, b1)
f1_1[n] = sinfunc(x1(n), a1, b1) #sin
f2_1[n] = expfunc(x1(n), a1, b1) #exp
f11_1[n] = linearfunc(x1(n), a1, b1) #linear
f21_1[n] = quadfunc(x1(n), a1, b1) #quad


for m in range(M_1 + 1):
E0_1[n] = g0(w1(m), a1, b1)
E1_1[n] = gsin(w1(m), a1, b1)
E2_1[n] = gexp(w1(m), a1, b1)
E11_1[n] = g1(w1(m), a1, b1)
E21_1[n] = g2(w1(m), a1, b1)

C_1[m, 0]= C0(w1(m), h_1, a)

C_1[m, M_1] = CN(w1(m), h_1, b)
for n in range(1,N1):
C_1[m, n] = Cb(n - 1, w1(m), h_1, a)

F0_1 = C_1 * np.transpose(f0_1)

F1_1 = C_1 * np.transpose(f1_1)
F2_1 = C_1 * np.transpose(f2_1)
F11_1 = C_1 * np.transpose(f11_1)
F21_1 = C_1 * np.transpose(f21_1)

# h = 0.001
w2 = np.arange(a,b+0.001,0.001)
M_2 = len(w2) - 1
h_2 = (b - a) / M_2
x2 = np.arange(a,b+0.001,0.001)
N2 = len(x2) - 1

C_2 = np.zeros(len(w2), len(x2))

E0_2 = np.zeros(1, len(w2))
E1_2 = np.zeros(1, len(w2))
E2_2 = np.zeros(1, len(w2))
E11_2 = np.zeros(1, len(w2))
E21_2 = np.zeros(1, len(w2))

f0_2 = np.zeros(1, len(x2))

f1_2 = np.zeros(1, len(x2))
f2_2 = np.zeros(1, len(x2))
f11_2 = np.zeros(1, len(x2))
f21_2 = np.zeros(1, len(x2))

for n in range(N2 + 1):

f0_2[n] = rectfunc(x2(n), a1, b1)
f1_2[n] = sinfunc(x2(n), a1, b1) # sin
f2_2[n] = expfunc(x2(n), a1, b1) # exp
f11_2[n] = linearfunc(x2(n), a1, b1) # linear
f21_2[n] = quadfunc(x2(n), a1, b1) # quad


for m in range(M_2 + 1):
E0_2[n] = g0(w2(m), a1, b1)
E1_2[n] = gsin(w2(m), a1, b1)
E2_2[n] = gexp(w2(m), a1, b1)
E11_2[n] = g1(w2(m), a1, b1)
E21_2[n] = g2(w2(m), a1, b1)

C_2[m, 0] = C0(w2(m), h_2, a)

C_2[m, M_2] = CN(w2(m), h_2, b)
for n in range(1,N2):
C_2[m, n] = Cb(n - 1, w2(m), h_2, a)

F0_2 = C_2 * np.transpose(f0_2)

F1_2 = C_2 * np.transpose(f1_2)
F2_2 = C_2 * np.transpose(f2_2)
F11_2 = C_2 * np.transpose(f11_2)
F21_2 = C_2 * np.transpose(f21_2)

###################
xx = np.arange(-2,2+0.001,0.001)
ff = np.zeros(1, len(xx))
fexp = np.zeros(1, len(xx))
fsin = np.zeros(1, len(xx))
flinear = np.zeros(1, len(xx))
fquad = np.zeros(1, len(xx))
for n in range(len(xx)):
ff[n] = rectfunc(xx(n), a1, b1)
fsin[n] = sinfunc(xx(n), a1, b1)
fexp[n] = expfunc(xx(n), a1, b1)
flinear[n] = linearfunc(xx(n), a1, b1)
fquad[n] = quadfunc(xx(n), a1, b1)

import matplotlib.pyplot as plt

fig, ax = plt.subplots(2, 2, sharex='col', sharey='row')
ax[0][0].plot(x,y1)
ax[0][1].plot(x,y2)
ax[1][0].plot(x,y3)
ax[1][1].plot(x,y4)
## OQF and FFT[f0]
# f11
subplot(3, 2, [1, 2])
plot(xx, ff)
axis([-1.5 1.5 - 0.5 1.5])
title('$f_0$', 'Interpreter', 'latex');

subplot(3, 2, 3)

plot(w, real(E0), w1, real(F0_1), w, real(F0