weight.m

用无网格法计算VAN DE POL 振子

function [w, dwdx, dwd2x, dwd3x] = Weight(type, para, di, dmi)
% EVALUATE WEIGHT FUNCTION
%
% SYNTAX: [w, dwdx, dwd2x, dwd3x, dwd4x] = Weight(type, para, di, dmi)
%
% INPUT PARAMETERS
%    type - Type of weight function
%    para - Weight function parameter
%    di   - Distance
%    dmi  - Support size
% OUTPUT PARAMETERS
%    w    - Value of weight function at r
%    dwdx - Value of first order derivative of weight function with respect to x at r
%    dwd2x- Value of Second order derivative of weight function with respect to x at r
%
r  = abs(di) / dmi;

if (di >= 0.0)
   drdx = 1.0/dmi;
else
   drdx = -1.0/dmi;
end

% EVALUATE WEIGHT FUNCTION AND ITS FIRST AND SECOND ORDER OF DERIVATIVES WITH RESPECT r AT r
if (type == 'GAUSS')
   [w,dwdr,dwd2r,dwd3r,dwd4r] = Gauss(para,r);
elseif (type == 'CUBIC')
   [w,dwdr,dwd2r,dwd3r,dwd4r] = Cubic(r);
elseif (type == 'SPLIN')
   [w,dwdr,dwd2r,dwd3r,dwd4r] = Spline(r);
elseif (type == 'CSRBF')
   [w,dwdr,dwd2r,dwd3r,dwd4r] = CSRBF2(r);
else
   error('Invalid type of weight function.');
end

w = w;
dwdx  = dwdr * drdx;
dwd2x = dwd2r * drdx.^2;
dwd3x = dwd3r * drdx.^3;
dwd4x = dwd4r * drdx.^4;


% SUBROUTINES
function [w,dwdr,dwd2r,dwd3r,dwd4r] = Gauss(beta,r)
if (r>1.0)
   w     = 0.0;
   dwdr  = 0.0;
   dwd2r = 0.0;
   dwd3r = 0.0;
   dwd4r = 0.0;
else
   b2 = beta*beta;
   r2 = r*r;
   eb2 = exp(-b2);

   w     = (exp(-b2*r2) - eb2) / (1.0 - eb2);
   dwdr  = -2*b2*r*exp(-b2*r2) / (1.0 - eb2);
   dwd2r = -2*b2*exp(-b2*r2)*(1-2*b2*r2) / (1.0 - eb2);
   dwd3r = 4*b2^2*r*exp(-r2*b2)*(-3+2*r2*b2)/(-1+exp(-b2)); 
   dwd4r = -4*b2^2*exp(-r2*b2)*(3-12*r2*b2+4*r2^2*b2^2)/(-1+exp(-b2));
end

function [w,dwdr,dwd2r,dwd3r,dwd4r] = Cubic(r)
if (r>1.0)
   w     = 0.0;
   dwdr  = 0.0;
   dwd2r = 0.0;
   dwd3r = 0.0;
   dwd4r = 0.0;
else
   w     = 1-6*r^2+8*r^3-3*r^4;
   dwdr  = -12*r+24*r^2-12*r^3;
   dwd2r = -12+48*r-36*r^2;
   dwd3r = 48-72*r;
   dwd4r = -72;
end

function [w,dwdr,dwd2r,dwd3r,dwd4r] = Spline(r)
if (r>1.0)
   w     = 0.0;
   dwdr  = 0.0;
   dwd2r = 0.0;
   dwd3r = 0.0;
   dwd4r = 0.0;
elseif (r<=0.5)
   w     = 2/3 - 4*r^2 + 4*r^3;
   dwdr  = -8*r + 12*r^2;
   dwd2r = -8 + 24*r;
   dwd3r = 24;
   dwd4r = 0;
else
   w     = 4/3 - 4*r + 4*r^2 - 4*r^3/3;
   dwdr  = -4 + 8*r -4*r^2;
   dwd2r = 8 - 8*r;
   dwd3r = -8;
   dwd4r = 0;
end

function [w,dwdr,dwd2r,dwd3r,dwd4r] = CSRBF2(r)
if (r>1.0)
   w     = 0.0;
   dwdr  = 0.0;
   dwd2r = 0.0;
   dwd3r = 0.0;
   dwd4r = 0.0;
else
   w  = (1-r)^6*(6+36*r+82*r^2+72*r^3+30*r^4+5*r^5);
   dwdr  = 11*r*(r+2)*(5*r^3+15*r^2+18*r+4)*(r-1)^5;
   dwd2r = 22*(25*r^5+100*r^4+142*r^3+68*r^2-16*r-4)*(r-1)^4;
   dwd3r = 4752*r+20790*r^4-13860*r^2-16632*r^6+4950*r^8;
   dwd4r = 4752-27720*r+83160*r^3-99792*r^5+39600*r^7;
end