Numerical solution of an integral equation 
March 2010, April 2016

We solve the problem of steady flow of a gas down a cylindrical tube of length L and radius R, which lies along the x axis from x=0 to L.  The gas is assumed to have such a low density that the molecules make collisions only with the walls of the tube and not with one another.  When they strike the wall, they stick to it for a short time and are then re-emitted in a random direction, so that the probability p(|x-y|) δ x of getting from point y on the wall to another point in the range [x,x+δx] can be calculated using only the geometry of the tube.  In a steady flow, the molecules must leave the surface at the same rate they arrive.  Accordingly, the rate ν(x) at which molecules leave point x on the surface satisfies the integral equation

({{({{ν(x)=s(x)+}})(∫_0^L)p(|x-y|)ν(y)dy}})

where s(x) is the rate at which molecules arrive at x from the open end of the tube at x=0 and the second term on the RHS is the rate of arrival of molecules that have been previously been absorbed and then re-emitted from the tube wall. 

First we find an expression for the probability (or rather, the probability density) for a molecule to travel a distance x between collisions with the walls of the tube.  Suppose the molecule starts from a point on the wall at x=0 and is emitted in a random direction, travelling a vector distance D to arrive at an element of surface dS at a distance x further down the tube.  (Note that D isn't the same as x, because although the molecule travels in a straight line, it needn't be travelling parallel to the x axis.)   You should be able to persuade yourself that, if the normal to the surface at the starting point is n, the probability of the molecule arriving at dS is

({{({{(((OverHat[D]) · n)/π) × (((OverHat[D]) · dS)/(D^2))}})}})

From the symmetry between the point of emission and the point of arrival, the angle between D and dS is the same as that between D and n, so the last expression can be simplified to

({{({{(((((OverHat[D]) · n))^2)/(π (D^2)))}})× dS}})

where dS = R d ϕ  dx and ϕ is the angle measured around the axis of the tube.  We just need to integrate with respect to ϕ to obtain p(x) dx.  

Let's obtain an explicit expression for the prefactor multiplying dS.  The molecule starts from

x0={0,0,R}

where the inward normal to the surface is

n={0,0,-1}

and it ends up at position

x1={x,R Sin[ϕ],R Cos[ϕ]}

So, calculating the remaining terms...

Dee=x1-x0

prefactor=FullSimplify[(Dee.n/Dee.Dee)^2]/π

From this point onward, all lengths will be measured in units of R, so we replace R by unity:

prefactor=prefactor/.R->1

This is the quantity that needs to be integrated with respect to ϕ to obtain p(x).  Although the angular integration is straightforward it is certainly messy and it's convenient to let Mathematica do the job. To make its work easier (and the output simpler), we specify that x is real and positive:

p = Integrate[prefactor,{ϕ,0,2π},Assumptions->{x>0}]

The probability density must integrate to unity, which we verify:

2Integrate[p,{x,0,Infinity}]

where the factor of 2 takes account of the integration from -∞ to 0.

Next we integrate p from x to ∞ to obtain the "source function" as a function of x , i.e., the rate of arrival of molecules from x<0.  The reasoning here uses the solution of the problem for the case L= ∞. In this limit, the net rate of flow down the tube is zero and the gas is at equilibrium throughout the system; hence the rate at which molecules leave the wall is independent of position, and equal to  (Subscript[ν, 0]), say: we choose our units so that (Subscript[ν, 0])=1. Substituting this value into the integral equation gives

({{({{1=s(x)+(∫_0^∞)p(|x-y|)dy}})}})

By using the symmetry of p(x) and the fact that it integrates to give 1,the last equation can be rearranged as s(x)=(∫_x^∞)p(y)dy. But since s(x) cannot depend on the length of the tube (remember that the molecules arrive from outside the tube), we can use this same expression in the general case.  Accordingly,

s=Integrate[p,{x,y,Infinity},Assumptions->{y>0}]/.y->x

Finally, we define a function flow to solve the integral equation

({{({{ν(x)=s(x)+}})(∫_0^L)p(|x-y|)ν(y)dy}})

and plot the emission rate ν(x) from the wall, as a function of position, for a given source function s(x) and transfer probability p(x). The method we use is to discretize x and y, approximating the integral by a simple sum of rectangles.  If we define (Subscript[x, i])=(i-0.5)L/n, (Subscript[s, i])=s((Subscript[x, i])),  (Subscript[ν, i])=ν((Subscript[x, i])) and (Subscript[K, i j])=L p(|(Subscript[x, i])-(Subscript[x, j])|)/n for i,j=1 to n, the integral equation becomes

({{({{(Subscript[ν, i])=(Subscript[s, i])+(Subscript[K, i j])(Subscript[ν, j])}})}})

or

({{({{(Subscript[M, i j])(Subscript[ν, j])=(Subscript[s, i])     where    (Subscript[M, i j])=(Subscript[δ, i j])-(Subscript[K, i j])}})}})

The input parameters to the function are the length of the tube (in units of R) and the number n of subdivisions used in discretizing the problem. The Mathematica function Module is used to define local variables {L,step,...} that are not needed elsewhere in the script; it's assumed that the expressions p and s have already been calculated, as above.

flow[LL_,n_]:=
Module[{L,step,xx,source,K,M,ν}, 
L =N[LL]; (* tube length in units of tube radius *)
step = L/n; (* width of one rectangle *)
xx=step Table[i-0.5,{i,n}]; (* x values are midpoints of intervals *)
source=Table[s/.x->xx[[i]],{i,n}]; (* source values in units of ν0 *)
K = step Table[p/.x->(Abs[xx[[i]]-xx[[j]]]),{i,n},{j,n}];(* kernel, represented as a matrix *)
M = IdentityMatrix[n]-K;
ν=LinearSolve[M,source]; (* emission rate in units of ν0 *)
ListPlot[{Table[{xx[[i]],source[[i]]},{i,n}],Table[{xx[[i]],ν[[i]]},{i,n}]},
PlotRange->{{0,L},{0,Full}},Joined->True,
PlotStyle->{{Blue,Dashed},Blue},
Frame->True,
FrameLabel->{{"ν/(Subscript[ν, 0]) and s/(Subscript[ν, 0])",None},{"x/R","emission rate and source (dashed) vs position"}},RotateLabel->True]]

As an example, here is the solution for a tube of length 5R, subdivided into 100 slices of length 0.05R :

flow[5,100]

If you want to change the input values you will first have to set up the expressions p and s, and the function flow, by selecting "Evaluate Notebook" from the "Evaluation" menu.  It may take a minute or two for Mathematica to do the integrations.

Note that the solution for the emission rate is remarkably close to the straight line ν(x)=(L+1-x)/(L+2), which would be the exact solution if the transfer probability had the model form p(x)=(1/2)(e^(-|x|)) (see Examples Sheet 3.9); the maximum deviation away from a straight line increases as L increases, but is never large.  You might like to modify this script  to investigate this behaviour and/or compare the numerical results with the exact result for the exponential p(x).