（个人学习笔记，仅供参考）

import numpy as np
from scipy.special import kv, erfc
from scipy.integrate import dblquad
import matplotlib.pyplot as plt
import scipy.integrate as spi

# Constants
w = 0.6198
g0_sq = 21.5989
rho = 0.782908
Q2 = 10.0
P2 = 0.01
z0 = 4.25

for x in x_values:
    result, error = compute_integral(x)
    integral_results178.append(result)
    errors.append(error)

# Convert results to numpy arrays for further processing or plotting
integral_results178 = np.array(integral_results178)
errors = np.array(errors)

integral_results178T = integral_results178 - integral_results178L

integral_results178R = integral_results178L/integral_results178T

plt.figure(figsize=(5, 3))
plt.plot(x_values, integral_results1R, label='$Q^2=1$GeV$^2$', color='red')
plt.plot(x_values, integral_results178R, label='$Q^2=10$GeV$^2$', color='blue')
plt.xlabel('x')
plt.xscale('log')
plt.xlim(1e-5, 0.1)
plt.ylim(-0.05, 0)  # Adjust y-axis based on max result
plt.ylabel('$F_2^{\gamma}(x,Q^2)$')
plt.legend()
plt.show()