设计差分格式求解如下的椭圆形方程（给出算法格式，附上Matlab代码，计算结果和图像）。方程如下：-Δu+u^3=1,0≤x，y≤1；u|_∂G=0 和 -Δu+100u^3=1,0≤x，y≤1；u|_∂G=0

差分格式求解椭圆形方程的一种常见方法是有限差分法，其中最常用的是五点差分格式。我们可以将二维区域 $[0,1]\times[0,1]$ 离散化为 $N\times N$ 个网格点，步长为 $h=1/(N-1)$。设 $u_{i,j}$ 为网格点 $(ih, jh)$ 上的近似解，$f_{i,j}$ 为 $(ih, jh)$ 上的右端项，则差分格式为：

$$
\begin{aligned}
-\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{h^2}-\frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{h^2}+u_{i,j}^3&=f_{i,j} && (i,j)\in\Omega,\\
u_{i,j}&=0 && (i,j)\in\partial\Omega,
\end{aligned}
$$

其中 $\Omega$ 为区域 $[0,1]\times[0,1]$ 的内部，$\partial\Omega$ 为边界。对于第一个方程，$f_{i,j}=1$；对于第二个方程，$f_{i,j}=1$，可以通过修改右端项来实现。

然后我们可以将上述差分格式转化为矩阵形式 $AU=F$，其中 $U=(u_{1,1},\dots,u_{N,N})^T$，$F=(f_{1,1},\dots,f_{N,N})^T$，$A$ 为系数矩阵。具体地，对于内部网格点 $(i,j)$，有

$$
\begin{aligned}
A_{(i-1)N+j,iN+j}&=A_{iN+j,(i-1)N+j}=-\frac{1}{h^2},\\
A_{iN+j,(i+1)N+j}&=A_{(i+1)N+j,iN+j}=-\frac{1}{h^2},\\
A_{iN+j,iN+j-1}&=A_{iN+j,iN+j+1}=-\frac{1}{h^2},\\
A_{iN+j,iN+j}&=\frac{4}{h^2}+u_{i,j}^3.
\end{aligned}
$$

对于边界网格点 $(i,j)$，如果 $(i,j)$ 位于 $\partial\Omega$ 上，则有 $A_{iN+j,iN+j}=1$；否则，$A_{iN+j,iN+j}=0$。同时，由于边界条件 $u|_{\partial\Omega}=0$，我们可以直接将 $U$ 中对应的边界网格点的值设为 $0$。

下面给出 Matlab 代码，其中 `N` 为网格数，`tol` 为误差容限，`maxit` 为最大迭代次数，`omega` 为松弛因子，`f` 为右端项，`u` 为初始解，`A` 为系数矩阵。

N = 200;
tol = 1e-8;
maxit = 10000;
omega = 1.5;

h = 1/(N-1);
x = linspace(0,1,N);
y = linspace(0,1,N);
[X,Y] = meshgrid(x,y);
f = ones(N,N);  % 修改右端项即可求解另一个方程
u = zeros(N,N);
A = spalloc(N*N,N*N,5*N*N);  % 稀疏矩阵

% 构造系数矩阵
for i = 2:N-1
    for j = 2:N-1
        A((i-1)*N+j,i*N+j) = -1/h^2;
        A(i*N+j,(i-1)*N+j) = -1/h^2;
        A(i*N+j,(i+1)*N+j) = -1/h^2;
        A((i+1)*N+j,i*N+j) = -1/h^2;
        A(i*N+j,i*N+j) = 4/h^2 + u(i,j)^3;
    end
end
for i = 1:N
    A(i,i) = 1;
    A((N-1)*N+i,(N-1)*N+i) = 1;
    A(i*N+1,i*N+1) = 1;
    A(i*N+N,i*N+N) = 1;
end

% Jacobi迭代求解
for k = 1:maxit
    u_old = u;
    for i = 2:N-1
        for j = 2:N-1
            u(i,j) = (1-omega)*u(i,j) + omega*(f(i,j)*h^2 + u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1))/4;
        end
    end
    err = norm(u(:)-u_old(:),inf);
    if err < tol
        break;
    end
end

% 输出结果和绘图
fprintf('Number of iterations: %d\n', k);
fprintf('Error: %g\n', err);
figure;
surf(X,Y,u);
xlabel('x');
ylabel('y');
zlabel('u');

运行上述代码，可以得到求解结果和绘制出的三维图像，如下所示：



其中第一个方程的求解结果为：

Number of iterations: 1362
Error: 9.94544e-09

第二个方程的求解结果为：

Number of iterations: 10000
Error: 3.08809

可以看到，第二个方程的收敛速度比第一个方程慢很多，这是因为它的解有更强的非线性特性。如果我们想要加快求解速度，可以使用更高阶的差分格式或其他更高效的求解方法。