给我3个谱元法求解二维非线性Klein-Gordon方程的matlab代码

以下是三种谱元法求解二维非线性Klein-Gordon方程的Matlab代码：

第一种方法：Legendre谱元法

clear all; clc;
N=31; M=31; c=1; a=1; dt=0.01; T=1;
x=linspace(-a,a,N); y=linspace(-a,a,M); [X,Y]=meshgrid(x,y);
u0=exp(-5*(X.^2+Y.^2)); v0=zeros(N,M);
for n=1:round(T/dt)
    u=u0+(dt/2)*v0;
    [L,U]=legDc(N-1,2*a); [Lp]=legDcDx(N-1,L);
    [P,W]=legDcWt(N-1,-a,a); W=diag(W);
    A=kron(eye(M),Lp)*kron(W,eye(N))+kron(U,W)*kron(eye(N),Lp);
    B=kron(eye(M),diag(1./P.^2))*kron(W,eye(N))+c^2*kron(diag(1./P.^2),W)*kron(eye(N),diag(1./P.^2));
    G=@(x) 2*sin(2*x); Gp=@(x) 4*cos(2*x);
    F=(kron(eye(M),diag(1./P.^2))*kron(W,eye(N))-c^2*kron(diag(1./P.^2),W)*kron(eye(N),diag(1./P.^2)))*u(:)-v0(:)+dt*Gp(u0(:)).*v0(:)-dt*G(u0(:));
    u=u+reshape(A\F(:),N,M);
    v=v0+(dt/2)*(Gp(u0).*v0-G(u0));
    u0=u; v0=v;
end
mesh(X,Y,u0);

第二种方法：Chebyshev谱元法

clear all; clc;
N=31; M=31; c=1; a=1; dt=0.01; T=1;
x=linspace(-a,a,N); y=linspace(-a,a,M); [X,Y]=meshgrid(x,y);
u0=exp(-5*(X.^2+Y.^2)); v0=zeros(N,M);
for n=1:round(T/dt)
    u=u0+(dt/2)*v0;
    [D,x]=cheb(N); D2=D^2; D2=D2(2:N,2:N);
    [P,w]=clencurt(N-1); W=diag(w);
    A=kron(eye(M),D2)*kron(W,eye(N))+kron(D2,W)*kron(eye(N),eye(M));
    B=kron(eye(M),diag(1./P.^2))*kron(W,eye(N))+c^2*kron(diag(1./P.^2),W)*kron(eye(N),diag(1./P.^2));
    G=@(x) 2*sin(2*x); Gp=@(x) 4*cos(2*x);
    F=(kron(eye(M),diag(1./P.^2))*kron(W,eye(N))-c^2*kron(diag(1./P.^2),W)*kron(eye(N),diag(1./P.^2)))*u(:)-v0(:)+dt*Gp(u0(:)).*v0(:)-dt*G(u0(:));
    u=u+reshape(A\F(:),N,M);
    v=v0+(dt/2)*(Gp(u0).*v0-G(u0));
    u0=u; v0=v;
end
mesh(X,Y,u0);

第三种方法：Fourier谱元法

clear all; clc;
N=31; M=31; c=1; a=1; dt=0.01; T=1;
x=linspace(-a,a,N); y=linspace(-a,a,M); [X,Y]=meshgrid(x,y);
u0=exp(-5*(X.^2+Y.^2)); v0=zeros(N,M);
for n=1:round(T/dt)
    u=u0+(dt/2)*v0;
    [X,Y]=meshgrid(x,y); kx=pi/a*[-floor(N/2):floor((N-1)/2)]; ky=pi/a*[-floor(M/2):floor((M-1)/2)];
    [KX,KY]=meshgrid(kx,ky); KX=KX(:); KY=KY(:);
    A=-c^2*(KX.^2+KY.^2); B=1./A; B(find(isnan(B)))=0;
    G=@(x) 2*sin(2*x); Gp=@(x) 4*cos(2*x);
    F=B.*u(:)-v0(:)+dt*Gp(u0(:)).*v0(:)-dt*G(u0(:));
    u=real(ifft2(reshape(F,N,M).*reshape(A,N,M)));
    v=v0+(dt/2)*(Gp(u0).*v0-G(u0));
    u0=u; v0=v;
end
mesh(X,Y,u0);