使用谱方法考虑求解如下的亥姆霍兹方程： $$\nabla^2 u(x,y) + k^2u(x,y) = f(x,y)$$，基函数选择傅里叶级数，给出Python程序

以下为使用谱方法求解二维亥姆霍兹方程的Python程序：

import numpy as np
import matplotlib.pyplot as plt

# 定义问题中的参数
Lx = 1.0  # 区域长度
Ly = 1.0  # 区域宽度
Nx = 64  # 空间离散点数
Ny = 64
k = 10.0  # 波数
f = lambda x, y: np.exp(-((x-0.5)**2 + (y-0.5)**2)/0.1)  # 源函数

# 计算傅里叶系数
def fourier_coeff(f, Nx, Ny):
    fx = np.zeros((Nx, Ny))
    fy = np.zeros((Nx, Ny))
    for j in range(Ny):
        fy[:, j] = np.fft.fft(f[:, j])
    for i in range(Nx):
        fx[i, :] = np.fft.fft(fx[i, :])
    return fx, fy

# 计算解函数
def solve_helmholtz(Lx, Ly, Nx, Ny, k, f):
    # 创建网格和求解矩阵
    x = np.linspace(0, Lx, Nx, endpoint=False)
    y = np.linspace(0, Ly, Ny, endpoint=False)
    X, Y = np.meshgrid(x, y, indexing='ij')
    dx = Lx / Nx
    dy = Ly / Ny
    A = np.zeros((Nx, Ny, Nx, Ny))
    for i in range(Nx):
        for j in range(Ny):
            A[i, j, i, j] = -4.0 + k**2 * dx**2
            if i > 0:
                A[i, j, i-1, j] = 1.0
            if i < Nx-1:
                A[i, j, i+1, j] = 1.0
            if j > 0:
                A[i, j, i, j-1] = 1.0
            if j < Ny-1:
                A[i, j, i, j+1] = 1.0
    A = A.reshape((Nx*Ny, Nx*Ny))
    
    # 计算傅里叶系数和右端项
    fx, fy = fourier_coeff(f(X, Y), Nx, Ny)
    rhs = np.zeros(Nx*Ny, dtype=np.complex)
    for i in range(Nx):
        for j in range(Ny):
            rhs[i*Ny+j] = fx[i, j] + fy[i, j]
    
    # 解线性方程组
    u = np.linalg.solve(A, rhs)
    u = u.reshape((Nx, Ny)).real
    
    return u

# 计算解函数
u = solve_helmholtz(Lx, Ly, Nx, Ny, k, f)

# 绘图
plt.imshow(u, extent=[0, Lx, 0, Ly], origin='lower')
plt.colorbar()
plt.show()