bkm2dcondk.m

％ A 2D homogeneous convection-diffusion case (u=exp(-ex*deta*x-ex*deta*y) with a square with ％ all

% A 2D homogeneous convection-diffusion case (u=exp(-ex*deta*x-ex*deta*y) with a square with 
% all Dirichlet boundary, note that reaction coefficient is not zero
% by indirect  BKM 

% Note that Sbesseli is a scaled Besseli function with scaler exp(-Lamda*a*c+Lamda*r)

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

   kd=deta^2*3/2;  % reaction coefficient
   Vx=-deta;       % velocity along x-axis
   Vy=-deta;       % velocity along y-axis
   Lamda=sqrt((Vx^2+Vy^2)/4+kd);      % Bessel function coefficient in general solution
   
 scale=Lamda*sqrt(1)*a;  % scale down the largest entry in interpolation matrix   
 
   ex=1/2+sqrt(deta^2+2*kd)/(2*deta);   % ex is modification coefficients due to reaction k
  
    ed=ex*deta;     % exact solution coefficient 
   
% generating the BKM interpolation matrix  

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

 for i=1:4*n-4
     bb(i,1)=exp(-ed*Rx(i)-ed*Ry(i));
     for j=1:4*n-4
         r=D2r(Rx(i),Rx(j),Ry(i),Ry(j));  %  Euclidean distance
         Lr=Lamda*r;
         detax=Vx*(Rx(i)-Rx(j))/2;
         detay=Vy*(Ry(i)-Ry(j))/2;
         Q(i,j)=Con2d(scale,0,Lr,detax,detay); % Dirchelet boundary entries
     end
 end
 
%maxQ=max(max(Q))
Matrix_order=4*n-4
Peclect_Number=a*sqrt(2)*sqrt(2)*ed
%Matirx_rank=rank(Q)
Con_number=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));
             Lr=Lamda*r;
             detax=Vx*(xt(i1)-Rx(k))/2;
             detay=Vy*(yt(j1)-Ry(k))/2;
             au=au+alpha(k)*Con2d(scale,0,Lr,detax,detay);
         end
         u=exp(-ed*xt(i1)-ed*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