Consider the 1-D expression for heat transfer in a bar defined by the relation²

     (1)

where the 1-D domain is bounded by 0 ≤ x ≤ L. The exact solution to this problem is

     (2)

with the exact derivative of the temperature given by

     (3)

In order to solve the 1-D problem, a multiquadric (MQ) radial basis function (RBF) is used 

     (4)

where r(x, x_j) is the radial (Euclidean) distance from the expansion point (x_j) to any point (x) , c is a shape parameter that controls the flatness of the RBF and is set by the user, and n is an integer. With n = 1, we retrieve the inverse multiquadric

     (5)

that will be used to solve Eq. (1). Other types of RBFs are available; the MQ is accurate and popular.

A global expansion for the 1-D temperature can be expressed as

     (6)

with the second derivative of the temperature given as

     (7)

Introducing the RBF expansion for the terms in the governing equation, and collocating at the interior points, we obtain

     (8)

At the boundaries, we collocate the RBF expansion to impose the boundary conditions

     (9)

Defining the operator

     (10)

we can now assemble into a fully populated matrix as,

     (11)

The solutions obtained using finite difference, finite volume, finite element, boundary element, and the meshless method are listed in Table 1 for 6 equally spaced nodes³ with T_o = 15 and T_L = 25, and L = 1. The interior nodes do not have to be uniformly spaced.

Table 1. Comparison of errors for interior temperatures i = 2,3,…N-1

The Maple code listing follows:

> restart:
   with(LinearAlgebra):with(plots):

# MESHLESS METHOD SOLUTION USING MULTIQUADRIC RADIAL BASIS
FUNCTIONS (RBF) il:=6:To:=15:TL:=25:L:=1:
>   x:=[0,1/5,2/5,3/5,4/5,1]:
>   S:=1000:n:=1:dx:=1/(il-1):
>   C:=Array(1..il,1..il):phi:=Array(1..il,1..il):d2phi:=Array(1..il,1..il):
b:=Vector(1..il):TM:=Vector(1..il):alpha:=Vector(1..il):
for i from 1 to il do
   for j from 1 to il do
      phi[i,j]:=(1+(x[i]-x[j])^2/(S*dx^2))^(n-3/2):
      d2phi[i,j]:=3*((x[j]-x[i])/20)^2/(4*((x[j]-x[i])^2/40+1)^(5/2) )-1/(40*((x[j]-x[i])^2/40+1)^(3/2)):
   end do:
end do:
>   for i from 2 to il-1 do    
>       for j from 1 to il do
>         C[i,j]:=d2phi[i,j]+phi[i,j];
            b[i]:=-x[i];
>         C[1,j]:=phi[1,j];
           C[il,j]:=phi[il,j];
         end do:
      end do:
      b[1]:=To:b[il]:=TL:
> #ConditionNumber(C);
>   alpha:=LinearSolve(convert(C,Matrix),b):
TM[1]:=To:TM[6]:=TL:
for i from 2 to il-1 do
   for j from 1 to il do
>         TM[i]:=TM[i]+alpha[j]*(1+(x[i]-x[j])^2/(S*dx^2))^(n-3/2);
 >     end do:
>   end do:
     evalf(TM);

                     

> TE:=To*cos(xx)+(TL+L-To*cos(L))/sin(L)*sin(xx)-xx:
> TE:=subs(TE):
> TE:=plot(TE,xx=0..1,color=blue,legend="Exact",thickness=3):
> MEM:=[seq([subs(x[i]),subs(TM[i])],i=1..6)]:
> T:=plots[pointplot](MEM,style=line,color=red,legend="MEM", thickness=3): 
  MEM:=plots[pointplot](MEM,color=red,legend="MEM",symbol=box, symbolsize=15):

>  plots[display](TE,MEM,T,axes=BOXED,title="Solution - MEM");

              

Additional examples for two-dimensional domains are described in the text, along with a chapter on the boundary element method. The meshless method is an interesting numerical approach that belongs to the family of weighted residual techniques. The matrix condition number is on the order of 10¹⁰ and can give surprisingly good results – however, the solution can fluctuate when repeatedly executed, eventually returning to the nearly correct solution; this is not an issue when using local assembly instead of the global assemble performed here⁴. The method can be used to form hybrid schemes, e.g., a finite element method can easily be linked with a meshless method to solve a secondary system of equations for problems involving large domains. Results are not sensitive to the location of the nodes; a random placement of points gives qualitatively similar results as a uniform placement.