以下代码对误差的计算方法是错的，应该是解析解与数值解的误差，请你检查并修正。clc
clear
close all

% 初始化参数
N = 64;
L = 1;
h = L/N;
x = 0:h:L;
y = 0:h:L;
[X, Y] = meshgrid(x, y);
phi_init = zeros(N+1, N+1);
f = 2 * pi^2 * sin(pi*X).*sin(pi*Y);
tol = 0.001;
max_iter = 10000;

% 使用多种迭代方法求解
method = {'Jacobi','Gauss Seidel','SOR','Multigrid'};
phi = cell(1, 4);
error_iter = cell(1, 4);
num_iters = zeros(1, 4);

for k = 1:length(method)
phi{k} = phi_init;
error_iter{k} = zeros(1, max_iter);

for cnt = 1:max_iter
switch method{k}
case 'Jacobi'
phi{k} = Jacobi(phi{k}, f, h);
case 'Gauss Seidel'
phi{k} = Gauss_Seidel(phi{k}, f, h);
case 'SOR'
phi{k} = SOR(phi{k}, f, h, 1.5);
case 'Multigrid'
phi{k} = V_Cycle(phi{k}, f, h);
end

r = residual(phi{k}, f, h);
error_iter{k}(cnt) = max(max(abs(r)));
if error_iter{k}(cnt) < tol
break;
end
end

num_iters(k) = cnt;
disp([method{k} '迭代次数: ' num2str(cnt)  '误差:' num2str(max(max(abs(r))) )]);

% 画出迭代过程中的误差变化
figure(k);
semilogy(1:cnt, error_iter{k}(1:cnt));
xlabel('迭代次数');
ylabel('误差');
title([method{k} '误差随迭代次数的变化']);
grid on;
end

% 画出所有方法迭代过程中的误差变化对比
figure;
colors ='rgbk';
for k = 1:length(method)
semilogy(1:num_iters(k), error_iter{k}(1:num_iters(k)), colors(k), 'DisplayName', method{k}); 
hold on;
end
xlabel('迭代次数');
ylabel('误差');
title('不同方法误差随迭代次数的变化对比');
legend('show');
grid on;

% 分别定义求解函数：雅克比迭代、高斯赛德尔迭代、松弛法迭代、多重网格
function res = Jacobi(phi, f, h)
N = length(phi) - 1;
res = phi;

for j = 2:N
for i = 2:N
res(i, j) = (1/4) * (phi(i+1, j) + phi(i-1, j) + phi(i, j+1) + phi(i, j-1) - h^2 * f(i, j));
end
end
end

function res = Gauss_Seidel(phi, f, h)
N = length(phi) - 1;
res = phi;

for j = 2:N
for i = 2:N
res(i, j) = (1/4) * (res(i-1, j) + phi(i+1, j) + phi(i, j+1) + res(i, j-1) -h^2 * f(i, j));
end
end
end

function res = SOR(phi, f, h, omega)
N = length(phi) - 1;
res = phi;

for j = 2:N
for i = 2:N
old = res(i, j);
res(i, j) = (1/4) * (phi(i+1, j) + res(i-1, j) + phi(i, j+1) + res(i, j-1) - h^2 * f(i, j));
res(i, j) = old * (1 - omega) + omega * res(i, j);
end
end
end

function phi = V_Cycle(phi, f, h)
phi = Gauss_Seidel(phi, f, h);
r = residual(phi, f, h);
rhs = restriction(r);
eps = zeros(size(rhs));

if length(eps) == 3
eps = smoothing(eps, rhs, 2 * h);
else
eps = V_Cycle(eps, rhs, 2 * h);
end

phi = phi + prolongation(eps);
phi = smoothing(phi, f, h);
end

function res = smoothing(phi, f, h)
N = length(phi) - 1;
res = phi;

for j = 2:N
for i = 2:N
res(i, j) = (phi(i+1, j) + res(i-1, j) + phi(i, j+1) + phi(i, j-1) - h^2 * f(i, j)) / 4;
end
end
end

function res = residual(phi, f, h)
N = length(phi) - 1;
res = zeros(N + 1);

for j = 2:N
for i = 2:N
res(i, j) = f(i, j) - (phi(i+1, j) - 2 * phi(i, j) + phi(i-1, j) + phi(i, j+1) - 2 * phi(i, j) + phi(i, j-1)) / h^2;
end
end

end

function res = restriction(r)
N = (length(r) - 1) / 2;
res = zeros(N + 1);

for j = 2:N
for i = 2:N
res(i, j) = (r(2*i-2, 2*j) + 2 * r(2*i-1, 2*j) + r(2*i, 2*j)) / 4;
end
end
end

function res = prolongation(eps)
N = (length(eps) - 1) * 2;
res = zeros(N + 1);

for j = 2:2:N
for i = 2:2:N
res(i, j) = (eps(i/2, j/2) + eps(i/2+1, j/2) + eps(i/2, j/2+1) + eps(i/2+1, j/2+1)) / 4;
end
end

for j = 1:2:N+1
for i = 1:2:N+1
res(i, j) = eps((i+1)/2, (j+1)/2);
end
end
for j = 2:2:N
for i = 1:2:N
res(i, j) = (eps((i+1)/2, j/2+1) + eps((i+1)/2, j/2)) / 2;
end
end

for j = 1:2:N
for i = 2:2:N
res(i, j) = (eps(i/2+1, (j+1)/2) + eps(i/2, (j+1)/2)) / 2;
end
end
end