求解微分方程-Uxx-Uyy+Ux+20yUy+U=f，我已经写出了代码，但误差太大，始终检查不出错误。
思路是：fish2d是用来求解-Uxx-Uyy=f类型的快速泊松方程求解器，因此方程化为-Uxx-Uyy=f-Ux-20yUy-U。用初值u0=0代入右端来迭代，每次求解时调用fish2d函数。真解已知，根据真解求出了右端项f。
以下为代码（其实是matlab）：
- 主函数
clc,clear,close all
a=0;
b=1;
c=0;
d=1;
N1=31;
N2=31;
h1=(b-a)/N1;
h2=(d-c)/N2;
x=a+[1:N1-1]*h1;
y=c+[1:N2-1]*h2;
[X,Y]=meshgrid(x,y);
%准备矩阵
e1=ones(N1-1,1);
e2=ones(N2-1,1);
e3=ones((N1-1)*(N2-1),1);
I1=spdiags(e1,0,N1-1,N1-1);
I2=spdiags(e2,0,N2-1,N2-1);
I=spdiags(e3,0,(N1-1)*(N2-1),(N1-1)*(N2-1));
K3=spdiags([-e1,e1],[-1,1],N1-1,N1-1);
K4=spdiags([-e2,e2],[-1,1],N2-1,N2-1);
%A、B、C、D用以生成系数矩阵
C=kron(K3,I2);
C=C/(2*h1); %u_x系数矩阵
D=kron(I1,K4);
D=D/(2*h2); %u_y系数矩阵
for i=1:N2-1
k=i;
while k<=(N1-1)*(N2-1)
D(:,k)=D(:,k)*20*y(i); %a2（x，y）=20y。将矩阵每一列*20*所在列数
k=k+(N2-1);
end
end
%右端项f,代入真解解出
f=zeros(N1-1,N2-1);
for i=1:N1-1
for j=1:N2-1
f(i,j)=20*y(j)*(10*x(i)*y(j)*exp(x(i)^(9/2))*(x(i)- 1) + 10*x(i)*exp(x(i)^(9/2))*(x(i)- 1)*(y(j) - 1)) - 20*x(i)*exp(x(i)^(9/2))*(x(i)- 1)...
- 20*y(j)*exp(x(i)^(9/2))*(y(j) - 1) + 10*x(i)*y(j)*exp(x(i)^(9/2))*(y(j)- 1) + 10*y(j)*exp(x(i)^(9/2))*(x(i)- 1)*(y(j)- 1)...
- 90*x(i)^(9/2)*y(j)*exp(x(i)^(9/2))*(y(j) - 1) + 10*x(i)*y(j)*exp(x(i)^(9/2))*(x(i)- 1)*(y(j) - 1) - (405*x(i)^8*y(j)*exp(x(i)^(9/2))...
*(x(i)-1)*(y(j)-1))/2-(495*x(i)^(7/2)*y(j)*exp(x(i)^(9/2))*(x(i)- 1)*(y(j)-1))/2+ 45*x(i)^(9/2)*y(j)*exp(x(i)^(9/2))*(x(i)- 1)*(y(j) - 1);
end
end
f=reshape(f',(N1-1)*(N2-1),1);%reshape按列取值
% 求解
u0=zeros((N1-1)*(N2-1),1);
k=1;
while k<=2
u1=fish2d(f-(C+D+I)*u0);
u0=u1;
k=k+1;
end
%画图
u1=reshape(u1,N2-1,N1-1);
mesh(X,Y,u1);
title('数值解');
figure
u_true=10*X.*Y.*(1-X).*(1-Y).*exp(X.^4.5);
mesh(X,Y,u_true);
title('真解');
figure
mesh(X,Y,u_true-u1);
title('误差');
2.fish2d函数
function u = fish2d(f)
% Poisson solver in 2D based on matlab fft
% square geometry, homogeneous Dirichlet boundary conditions
% requires uniform mesh (2^k-1)^2 interior nodes
%
% C. T. Kelley, January 6, 1994
%
% This code comes with no guarantee or warranty of any kind.
%
% Solves - u_xx - u_yy = f
%
% Requires sintv.m and isintv.m
%
% function u = fish2d(f)
%
%
global fish_lal fish_ual ua;
n=length(f);
m=sqrt(n);
hh=.5/(m+1);
h2=(m+1)*(m+1);
e=ones(m,1);
%
% Don't compute eigenvalues, factor tridiagonal matrices
% or allocate space for the 2-D representation of the solution
% more than once. Note the use of global variables to store
% fish_ual, fish_lal, and ua.
%
if exist('fish_ual')==1
[f1, f2]=size(fish_ual);
m1=sqrt(f1);
else
m1=m;
end
if (exist('fish_ual') == 0) | (m1 ~= m)
%
% allocate room for ua, the 2D representation of the solution
%
ua=zeros(m,m);
%
% set up diagonal of system using eigenvalues of 1D operator
%
x=1:m; x=pi*hh*x'; lam=sin(x); lam=h2*(2*ones(m,1) + 4*lam.*lam);
en=h2*ones(n,1); fn=en; gn=en; dx=zeros(n,1);
for i=1:m
fn(i*m)=0;
gn((i-1)*m+1)=0;
dx((i-1)*m+1:i*m)=lam(i)*ones(m,1);
end
al=spdiags([-fn, dx, -gn], [-1,0,1], n,n);
%
% factor the tridiagonal matrix for solution in the x direction
%
[fish_lal,fish_ual]=lu(al);
end
%
% Map f into a rectangular array.
% Take sine transform in y direction.
%
ua(:)=f;
ua=sintv(ua).';
%
% Solve in the x direction
%
ua(:)=fish_ual\(fish_lal\ua(:));
%
% Take inverse transform.
%
ua=isintv((ua.'));
%
% Map uz into a linear array.
%
u=ua(:);
function u=sintv(z)
% sine transform, assume that size(z)=[nx,ny]=[2^k-1,ny]
%
% C. T. Kelley, May 1994
%
% This code comes with no guarantee or warranty of any kind.
%
[nx,ny]=size(z);
ww=-4*fft([zeros(ny,1), z']',2*nx+2);
u=imag(ww(2:nx+1,:));
function u=isintv(z)
% inverse sine transform, assume that size(z)=[nx,ny]=[2^k-1,ny]
%
% C. T. Kelley, May 1994
%
% This code comes with no guarantee or warranty of any kind.
%
[nx,ny]=size(z);
ww=ifft([zeros(ny,1), z']',2*nx+2);
u=imag(ww(2:nx+1,:));
