sbkm2d.m

A 2D homogeneous Helmholtz case (u=sin(x)cos(y) with a square) with ％ two Dirichlet edges (x=1,y=1

% A 2D homogeneous Helmholtz case (u=sin(x)cos(y) with a square) with 
% two Dirichlet edges (x=1,y=1) and two Neumann edges (x=0,y=0)
% by indirect symmetric BKM

  clear

  
  disp('Homogeneous 2D Helmholtz case by indirect symmetric BKM')
  
   n=input('input node number of one direction  '); %the node number of each direction
  
  t=clock;  % set initial time
  
  a=8; % length of rectangular
  b=8; % width of rectqangular
     
  [Rx,Ry]=Bsquare(a,b,n); % uniform boundary node generation
  
% Lamda value of Helmholtz equation

  Lamda=sqrt(2);
  
% generating the BKM interpolation matrix  

  Q=zeros(4*n-4,4*n-4); % preallocating memory
  bb=zeros(4*n-4,1);  

 % Dirichlet boundary 
 for i=1:2*n-2
     bb(i,1)=sin(Rx(i))*cos(Ry(i));
     for j=1:4*n-4
         r=D2r(Rx(i),Rx(j),Ry(i),Ry(j));  %  Euclidean distance
         if(j<=2*n-2)
         Q(i,j)=besselj(0,Lamda*r); % Dirchelet boundary entries
     elseif(j<=3*n-3)
         Q(i,j)=besselj(1,Lamda*r)*Lamda*(Rx(i)-Rx(j))/r; % Neumann entries at x=0
     else
         Q(i,j)=besselj(1,Lamda*r)*Lamda*(Ry(i)-Ry(j))/r; % Neumann entries at y=0
     end
 end
end


 % Neumann boundary at x=0 
 for i=2*n-2+1:3*n-3
     bb(i,1)=cos(Ry(i));
     for j=1:4*n-4
         r=D2r(Rx(i),Rx(j),Ry(i),Ry(j));  %  Euclidean distance
         if(j<=2*n-2)
         Q(i,j)=-besselj(1,Lamda*r)*Lamda*(Rx(i)-Rx(j))/r; % Dirchelet boundary entries
     elseif(j<=3*n-3&j==i)
         Q(i,j)=0.5*Lamda^2;
     elseif(j<=3*n-3&j~=i)
         Q(i,j)=-(-besselj(0,Lamda*r)+besselj(1,Lamda*r)/(Lamda*r))*(Rx(i)-Rx(j))^2 ...
             *Lamda^2/r^2+besselj(1,Lamda*r)*(Ry(i)-Ry(j))^2*Lamda/r^3; 
         % Neumann entries at x=0
     else
         Q(i,j)=-(-besselj(0,Lamda*r)+2*besselj(1,Lamda*r)/(Lamda*r))*(Rx(i)-Rx(j))...
             *(Ry(i)-Ry(j))*Lamda^2/r^2; 
         % Neumann entries at y=0
     end
 end
end


 % Neumann boundary at y=0 
 for i=3*n-3+1:4*n-4
     bb(i,1)=0;
     for j=1:4*n-4
         r=D2r(Rx(i),Rx(j),Ry(i),Ry(j));  %  Euclidean distance
         if(j<=2*n-2)
         Q(i,j)=-besselj(1,Lamda*r)*Lamda*(Ry(i)-Ry(j))/r; % Dirchelet boundary entries
     elseif(j<=3*n-3)
         Q(i,j)=-(-besselj(0,Lamda*r)+2*besselj(1,Lamda*r)/(Lamda*r))*(Ry(i)-Ry(j))...
             *(Rx(i)-Rx(j))*Lamda^2/r^2; 
         % Neumann entries at x=0
     elseif(j>3*n-3&j==i)
         Q(i,j)=0.5*Lamda^2;
     else
         Q(i,j)=-(-besselj(0,Lamda*r)+besselj(1,Lamda*r)/(Lamda*r))*(Ry(i)-Ry(j))^2 ...
             *Lamda^2/r^2+besselj(1,Lamda*r)*(Rx(i)-Rx(j))^2*Lamda/r^3; 
         % Neumann entries at y=0         
     end
 end
end

4*n-4
rank(Q)
condnumer=cond(Q)
alpha=Q\bb; % solving the RBF equation
clear bb;

% evaluating some BKM solutions in boundary and domain

  nt=10;  % nt^2 is the node number of interests
  [xt,yt]=square_eq(a,b,nt); % node co-ordinates of interests
  AerL2=0; % absolute error
  RerL2=0;  % relative error
  Aerr=size(nt,nt); %preallocating memory
  Rerr=size(nt,nt);
 for i1=1:nt
     for j1=1:nt
         au=0;
         for k=1:4*n-4
             r=D2r(xt(i1),Rx(k),yt(j1),Ry(k));
         if(k<=2*n-2)
             au=au+alpha(k)*besselj(0,Lamda*r); % Dirchelet boundary entries
         elseif(k<=3*n-3&r>=10^-6)
             au=au+alpha(k)*besselj(1,Lamda*r)*Lamda*(xt(i1)-Rx(k))/r; 
             % Neumann entries at x=0
         elseif(k<=4*n-4&r>=10^-6)
             au=au+alpha(k)*besselj(1,Lamda*r)*Lamda*(yt(j1)-Ry(k))/r; 
             % Neumann entries at y=0
         end
     end
     u=sin(xt(i1))*cos(yt(j1));  % Analytical solutions
     Aerr(i1,j1)=u-au;  % absolute error
     if(abs(u)<10^-3)  % relative error for aboslute value <0.001
         Rerr(i1,j1)=au;  % relative error
     else
         Rerr(i1,j1)=(u-au)/u;  % relative error
     end
     AerL2=AerL2+Aerr(i1,j1)^2; % square sum of absolute errors 
     RerL2=RerL2+Rerr(i1,j1)^2;   % square sum of relative errors
 end
end

AerL2=sqrt(AerL2/nt^2)   % L2 absolute error
RerL2=sqrt(RerL2/nt^2)   % L2 relative error

etime(clock,t)  % estimate the elapsed computing time (seconds)

 % err_helmholtz(n,xt,yt,Aerr,Rerr);  % show error disribution in 3D figure