7.3.1 plate_f1.m

% plate_f1.m -- Eigenvalue problem of square plates with free edges
%  written by Hirano Yuuki, Meiji University, Feb 2012
%  comments are modified by Masashi Katsurada, 24 June 2012.
%
%  \triangle^2 u=\lambda u
%  0<x<1, 0<y<1
%  mu: Poisson's ratio
%  usage:
%    A=plate_f1(640,0.3);
%    [v,d]=eigs(A,200,0);
%    plot_n(v(:,197),640,640)

function P1=plate_f1(N,mu)
        h=1/N;
        n=N+1;
        a=-2*(mu*mu+2*mu-3);
        b=1-mu*mu;
        c=-2*(mu-1);
        d=15-8*mu-5*mu*mu;
        e=-4*(mu*mu+mu-2);
        f=2-mu;
        g=-2*(mu-3);
        k=-2*(3*mu*mu+4*mu-8);
        In=speye(n,n);
        Jn=sparse(diag(ones(n-1,1),1)+diag(ones(n-1,1),-1));
        J2n=sparse(diag(ones(n-2,1),2)+diag(ones(n-2,1),-2));
        j2n=sparse(diag(ones(n-2,1),2)+diag(ones(n-2,1),-2));
        j2n(1,3)=sqrt(2);
        j2n(3,1)=sqrt(2);
        j2n(n-2,n)=sqrt(2);
        j2n(n,n-2)=sqrt(2);
        An=In;
        An(1,1)=b;
        An(n,n)=b;
        Bn=-8*In+2*Jn;
        Bn(1,1)=-e;
        Bn(n,n)=-e;
        Bn(1,2)=sqrt(2)*f;
        Bn(2,1)=sqrt(2)*f;
        Bn(n-1,n)=sqrt(2)*f;
        Bn(n,n-1)=sqrt(2)*f;
        Cn=-g*In+f*Jn;
        Cn(1,1)=-a;
        Cn(n,n)=-a;
        Cn(1,2)=sqrt(2)*f;
        Cn(2,1)=sqrt(2)*c;
        Cn(n-1,n)=sqrt(2)*c;
        Cn(n,n-1)=sqrt(2)*f;
        Dn=20*In-8*Jn+j2n;
        Dn(1,1)=k;
        Dn(n,n)=k;
        Dn(2,2)=19;
        Dn(n-1,n-1)=19;
        Dn(1,2)=-sqrt(2)*g;
        Dn(2,1)=-sqrt(2)*g;
        Dn(n-1,n)=-sqrt(2)*g;
        Dn(n,n-1)=-sqrt(2)*g;
        DDn=19*In-8*Jn+j2n;
        DDn(1,1)=d;
        DDn(n,n)=d;
        DDn(2,2)=18;
        DDn(n-1,n-1)=18;
        DDn(1,2)=-sqrt(2)*g;
        DDn(2,1)=-sqrt(2)*g;
        DDn(n-1,n)=-sqrt(2)*g;
        DDn(n,n-1)=-sqrt(2)*g;
        En=k*In-e*Jn+b*j2n;
        En(1,1)=2*a;
        En(n,n)=2*a;
        En(2,2)=d;
        En(n-1,n-1)=d;
        En(1,2)=-sqrt(2)*a;
        En(2,1)=-sqrt(2)*a;
        En(n-1,n)=-sqrt(2)*a;
        En(n,n-1)=-sqrt(2)*a;
        P1=kron(j2n,An)+kron(Jn,Bn)+kron(In,Dn);
        P1(1:n,1:n)=En;
        P1(1:n,n+1:2*n)=sqrt(2)*Cn';
        P1(n+1:2*n,1:n)=sqrt(2)*Cn;
        P1(n+1:2*n,n+1:2*n)=DDn;
        P1(n*(n-2)+1:n*(n-1),n*(n-2)+1:n*(n-1))=DDn;
        P1(n*(n-2)+1:n*(n-1),n*(n-1)+1:n*n)=sqrt(2)*Cn;
        P1(n*(n-1)+1:n*n,n*(n-2)+1:n*(n-1))=sqrt(2)*Cn';
        P1(n*(n-1)+1:n*n,n*(n-1)+1:n*n)=En;
        P1=P1/(h*h*h*h);
end