Note on the Solution of Transport Equation by Tau Method and Walsh Functions

Abstract

We consider the combined Walsh function for the three-dimensional case. A method for the solution of the neutron transport equation in three-dimensional case by using the Walsh function, Chebyshev polynomials, and the Legendre polynomials are considered. We also present Tau method, and it was proved that it is a good approximate to exact solutions. This method is based on expansion of the angular flux in a truncated series of Walsh function in the angular variable. The main characteristic of this technique is that it reduces the problems to those of solving a system of algebraic equations; thus, it is greatly simplifying the problem.

1. Introduction

The Walsh functions have many properties similar to those of the trigonometric functions. For example, they form a complete, total collection of functions with respect to the space of square Lebesgue integrable functions. However, they are simpler in structure to the trigonometric functions because they take only the values 1 and −1. They may be expressed as linear combinations of the Haar functions [1], so many proofs about the Haar functions carry over to the Walsh system easily. Moreover, the Walsh functions are Haar wavelet packets. For a good account of the properties of the Haar wavelets and other wavelets, see [2]. We use the ordering of the Walsh functions due to Paley [3, 4]. Any function 𝑓∈𝐿²[0,1] can be expanded as a series of Walsh functions𝑓(𝑥)=^∞_𝑖=0𝑐_𝑖𝑊_𝑖(𝑥),where𝑐_𝑖=¹₀𝑓(𝑥)𝑊_𝑖(𝑥). (1.1) In [5], Fine discovered an important property of the Walsh Fourier series: the 𝑚=2^𝑛th partial sum of the Walsh series of a function 𝑓 is piecewise constant, equal to the 𝐿¹ mean of 𝑓, on each subinterval [(𝑖−1)/𝑚,𝑖/𝑚]. For this reason, Walsh series in applications are always truncated to 𝑚=2^𝑛 terms. In this case, the coefficients 𝑐_𝑖 of the Walsh (-Fourier) series are given by 𝑐_𝑖=^𝑚−1_𝑗=01𝑚𝑊_𝑖𝑗𝑓_𝑗, (1.2) where 𝑓_𝑗 is the average value of the function 𝑓(𝑥) in the 𝑗th interval of width 1/𝑚 in the interval (0,1), and 𝑊_𝑖𝑗 is the value of the 𝑖th Walsh function in the 𝑗th subinterval. The order 𝑚 Walsh matrix, 𝒲_𝑚, has elements 𝑊_𝑖𝑗.

Let 𝑓(𝑥) have a Walsh series with coefficients 𝑐_𝑖 and its integral from 0 to 𝑥 have a Walsh series with coefficients 𝑏_𝑖: ∫^𝑥₀𝑓(𝑡)𝑑𝑡=∑^∞_𝑖=0𝑏_𝑖𝑊_𝑖(𝑥). If we truncate to 𝑚=2^𝑛 terms and use the obvious vector notation, then integration is performed by matrix multiplication 𝐛=𝑃^𝑇_𝑚𝐜, where𝑃^𝑇_𝑚=𝑃_𝑚/21𝐼2𝑚^𝑚/2−1𝐼2𝑚^𝑚/2𝑂_𝑚/2,𝑃^𝑇₂=1214−140, (1.3) and 𝐼_𝑚 is the unit matrix, 𝑂_𝑚 is the zero matrix (of order 𝑚), see [6].

2. The Three-Dimensional Spectral Solution

In the literature there several works on driving a suitable model for the transport equation in 2 and 3-dimensional case as well as in cylindrical domain, for example, see [7], and by using the eigenvalue error estimates for two-dimensional neutron transport, see [8], by applying the finite element method in an infinite cylindrical domain, see [9], similarly by using Chebyshev spectral-𝑆_𝑁 method, see [10], and the discrete ordinates in the infinite cylindrical domain, see [11].

In this paper, we consider combined Walsh function with the Sumudu transform in order to extend the transport problem for the three-dimensional case by following the similar method that was proposed in [7]. This method is based on expansion of the angular flux in a truncated series of Walsh function in the angular variable. By replacing this development in the transport equation, this will result a first-order linear differential system. First of all we consider the three-dimensional linear, steady state, transport equation which is given by 𝜇𝜕Ψ+𝜕𝑥𝐱,𝜇,𝜃1−𝜇²𝜕cos𝜃Ψ𝜕𝜕𝑦𝐱,𝜇,𝜃+sin𝜃Ψ𝜕𝑧𝐱,𝜇,𝜃+𝜎_𝑡Ψ=𝐱,𝜇,𝜃_1−102𝜋𝜎_𝑠𝜇^′,𝜃^′⟶𝜇,𝜃Ψ𝐱,𝜇^′,𝜃^′𝑑𝜃^′𝑑𝜇^′+𝑆𝐱,𝜇,𝜃, (2.1) where we assume that the spatial variable 𝐱∶=(𝑥,𝑦,𝑧) varies in the cubic domain Ω∶={(𝑥,𝑦,𝑧)∶−1≤𝑥,𝑦,𝑧≤1}, and Ψ(𝐱,𝜇,𝜃)∶=Ψ(𝑥,𝑦,𝑧,𝜇,𝜃) is the angular flux in the direction defined by 𝜇∈[−1,1] and 𝜃∈[0,2𝜋]. 𝜎_𝑡 and 𝜎_s denote the total and the differential cross section, respectively, 𝜎_𝑠(𝜇^′,𝜙^′→𝜇,𝜙) describes the scattering from an assumed pre-collision angular coordinates (𝜇^′,𝜃^′) to a postcollision coordinates (𝜇,𝜃) and 𝑆 is the source term. See [12] for further details.

Note that, in the case of one-speed neutron transport equation; taking the angular variable in a disc, this problem will corresponds to a three dimensional case with all functions being constant in the azimuthal direction of the 𝑧 variable. In this way the actual spatial domain may be assumed to be a cylinder with the cross-section Ω and the axial symmetry in 𝑧. Then 𝐷 will correspond to the projection of the points on the unit sphere (the “speed”) onto the unit disc (which coincides with 𝐷). See [13] for the details.

Given the functions 𝑓₁(𝑦,𝑧,𝜇,𝜙), 𝑓₂(𝑥,𝑧,𝜇,𝜙), and 𝑓₃(𝑥,𝑦,𝜇,𝜙) describing the incident flux, we seek for a solution of (2.1) subject to the following boundary conditions.

For the boundary terms in 𝑥, for 0≤𝜃≤2𝜋, letΨ=𝑓𝑥=±1,𝑦,𝑧,𝜇,𝜃¹𝑦,𝑧,𝜇,𝜃,𝑥=−1,0<𝜇≤1,0,𝑥=1,−1≤𝜇<0. (2.2) For the boundary terms in 𝑦 and for −1≤𝜇<1,Ψ=𝑓𝑥,𝑦=±1,𝑧,𝜇,𝜃²𝑥,𝑧,𝜇,𝜃,𝑦=−1,0<cos𝜃≤1,0,𝑦=1,−1≤cos𝜃<0. (2.3) Finally, for the boundary terms in 𝑧, for −1≤𝜇<1,Ψ=𝑓𝑥,𝑦,𝑧=±1,𝜇,𝜃³𝑥,𝑦,𝜇,𝜃,𝑧=−1,0≤𝜃<𝜋,0,𝑧=1,𝜋<𝜃≤2𝜋. (2.4)

Theorem 2.1. Consider the integrodifferential equation (2.1) under the boundary conditions (2.2), (2.3) and (2.4), then the function Ψ(𝑥,𝑦,𝑧,𝜇,𝜃) satisfy the following first-order linear differential equation system for the spatial component Ψ_𝑖,𝑗(𝑥,𝜇,𝜃)𝜇𝜕Ψ_𝑖,𝑗𝜕𝑥𝑥,𝜇,𝜃+𝜎_𝑡Ψ_𝑖,𝑗𝑥,𝜇,𝜃=𝐺_𝑖,𝑗𝑥;𝜇,𝜃¹₋₁𝜎_𝑠𝜇^′,𝜃^′⟶𝜇,𝜃Ψ_𝑖,𝑗𝑥,𝜇^′,𝜃^′𝑑𝜃^′𝑑𝜇^′, (2.5) with the boundary conditions Ψ_𝑖,𝑗−1,𝜇,𝜂=𝑓₁^𝑖,𝑗𝜇,𝜃, (2.6) where 𝑓₁^𝑖,𝑗=4𝜇,𝜃𝜋²¹₋₁𝑇_𝑖𝑦𝑅_𝑗(𝑧)1−𝑦²1−𝑧²𝑓₁Ψ𝑦,𝑧,𝜇,𝜃𝑑𝑧𝑑𝑦,_𝑖,𝑗𝐺1,−𝜇,𝜃=0,_𝑖,𝑗𝑥;𝜇,𝜃=𝑆_𝑖,𝑗𝑥,𝜇,𝜃−1−𝜇²cos𝜃^𝐼_𝑘=𝑖+1𝐴^𝑘_𝑖Ψ_𝑘,𝑗𝑥,𝜇,𝜃+sin𝜃^𝐽_𝑙=𝑗+1𝐵^𝑙_𝑗Ψ_𝑖,𝑙𝑥,𝜇,𝜃, (2.7) with 𝑆_𝑖,𝑗=4𝑥,𝜇,𝜃𝜋²¹₋₁𝑇_𝑖𝑦𝑅_𝑗(𝑧)1−𝑦²1−𝑧²𝑆𝐴𝐱,𝜇,𝜃𝑑𝑧𝑑𝑦,_𝑘𝑖=2𝜋¹₋₁𝑑𝑇𝑑𝑦^𝑘𝑦𝑇_𝑖𝑦1−𝑦²𝐵𝑑𝑦,_𝑙𝑗=2𝜋¹₋₁𝑑𝑅𝑑𝑦^𝑙𝑦𝑅_𝑗(𝑧)√1−𝑧²𝑑𝑧. (2.8)

Proof. Expanding the angular flux Ψ(𝑥,𝑦,𝑧,𝜇,𝜙) in a truncated series of Chebyshev polynomials 𝑇_𝑖(𝑦) and 𝑅_𝑗(𝑧) leads to Ψ𝑥,𝑦,𝑧,𝜇,𝜃=^𝐼𝐽_𝑖=0𝑗=0Ψ_𝑖,𝑗𝑥,𝜇,𝜃𝑇_𝑖𝑦𝑅_𝑗(𝑧). (2.9) We insert Ψ(𝑥,𝑦,𝑧,𝜇,𝜃) given by (2.9) into the boundary condition in (2.3), for 𝑦=±1. Multiplying the resulting expressions by 𝑅_𝑗(𝑧)/√1−𝑧² and integrating over 𝑧, we get the components Ψ_0,𝑗(𝑥,𝜇,𝜃) for 𝑗=0,…,𝐽, Ψ_0,𝑗𝑥,𝜇,𝜃=𝑓^𝑗₂𝑥,𝜇,𝜃−^𝐼_𝑖=1(−1)^𝑗Ψ_𝑖,𝑗Ψ𝑥,𝜇,𝜃,0<cos𝜃≤1,_0,𝑗𝑥,𝜇,𝜃=−^𝐼_𝑖=1Ψ_𝑖,𝑗𝑥,𝜇,𝜃,−1≤cos𝜃<0. (2.10) Similarly, we substitute Ψ(𝑥,𝑦,𝑧,𝜇,𝜃) from (2.9) into the boundary conditions for 𝑧=±1, multiply the resulting expression by 𝑇_𝑖(𝑦)/1−𝑦², 𝑖=0,…,𝐼 and integrating over 𝑦, to define the components Ψ_𝑖,0(𝑥,𝜇,𝜃). For −1≤𝑥≤1, −1<𝜇<1, Ψ_𝑖,0𝑥,𝜇,𝜃=𝑓^𝑖₃𝑥,𝜇,𝜃−^𝐽_𝑗=1(−1)^𝑗Ψ_𝑖,𝑗Ψ𝑥,𝜇,𝜃,0≤𝜃<𝜋,_𝑖,0𝑥,𝜇,𝜃=−^𝐽_𝑗=1Ψ_𝑖,𝑗𝑥,𝜇,𝜃,𝜋<𝜃≤2𝜋, (2.11) where 𝑓^𝛽₂𝑥,𝜇,𝜃=2−𝛿_0,𝑗𝜋¹₋₁𝑓₂𝑅𝑥,𝑧,𝜇,𝜃^𝑗(𝑧)√1−𝑧²𝑓𝑑𝑧,_𝑖3𝑥,𝜇,𝜃=2−𝛿_𝑖,0𝜋¹₋₁𝑓₃𝑇𝑥,𝑦,𝜇,𝜃^𝑖𝑦1−𝑦²𝑑𝑦. (2.12) To determine the components Ψ_𝑖,𝑗(𝑥,𝜇,𝜃), 𝑖=1,…,𝐼, and 𝑗=1,…,𝐽 we substitute Ψ(𝑥,𝜇,𝜃), from (2.3) into (2.1) and the boundary conditions for 𝑥=±1. Multiplying the resulting expressions by (𝑇_𝑖(𝑦)/1−𝑦²)×(𝑅_𝑗(𝑧)/√1−𝑧²), and integrating over 𝑦 and 𝑧 we obtain 𝐼×𝐽 one-dimensional transport problems, namely, 𝜇𝜕Ψ_𝑖,𝑗𝜕𝑥𝑥,𝜇,𝜃+𝜎_𝑡Ψ_𝑖,𝑗𝑥,𝜇,𝜃=𝐺_𝑖,𝑗𝑥;𝜇,𝜃¹₋₁𝜎_𝑠𝜇^′,𝜃^′⟶𝜇,𝜃Ψ_𝑖,𝑗𝑥,𝜇^′,𝜃^′𝑑𝜃^′𝑑𝜇^′, (2.13) with the boundary conditions Ψ_𝑖,𝑗−1,𝜇,𝜂=𝑓₁^𝑖,𝑗𝜇,𝜃, (2.14) where 𝑓₁^𝑖,𝑗=4𝜇,𝜃𝜋²¹₋₁𝑇_𝑖𝑦𝑅_𝑗(𝑧)1−𝑦²1−𝑧²𝑓₁Ψ𝑦,z,𝜇,𝜃𝑑𝑧𝑑𝑦,_𝑖,𝑗1,−𝜇,𝜃=0, (2.15) for 0<𝜇≤1, and 0≤𝜃≤2𝜋. Finally, 𝐺_𝑖,𝑗𝑥;𝜇,𝜃=𝑆_𝑖,𝑗𝑥,𝜇,𝜃−1−𝜇²cos𝜃^𝐼_𝑘=𝑖+1𝐴^𝑘_𝑖Ψ_𝑘,𝑗𝑥,𝜇,𝜃+sin𝜃^𝐽_𝑙=𝑗+1𝐵^𝑙_𝑗Ψ_𝑖,𝑙𝑥,𝜇,𝜃, (2.16) with 𝑆_𝑖,𝑗=4𝑥,𝜇,𝜃𝜋²¹₋₁𝑇_𝑖𝑦𝑅_𝑗(𝑧)1−𝑦²1−𝑧²𝑆𝐴𝐱,𝜇,𝜃𝑑𝑧𝑑𝑦,_𝑘𝑖=2𝜋¹₋₁𝑑𝑇𝑑𝑦^𝑘𝑦𝑇_𝑖𝑦1−𝑦²𝐵𝑑𝑦,_𝑙𝑗=2𝜋¹₋₁𝑑𝑅𝑑𝑦^𝑙𝑦𝑅_𝑗(𝑧)√1−𝑧²𝑑𝑧. (2.17) Now, starting from the solution of the problem given by (2.13)–(2.17) for Ψ_𝐼,𝐽(𝑥,𝜇,𝜃), we then solve the problems for the other components, in the decreasing order in 𝑖 and 𝑗. Recall that ∑^𝐼_𝑖=𝐼+1⋯=∑^𝐽_𝑗=𝐽+1≡0. Hence, solving 𝐼×𝐽 one-dimensional problems, the angular flux Ψ(𝐱,𝜇,𝜃) is now completely determined through (2.9).

Remark 2.2. If we deal with different type of boundary conditions, then we consider the first components Ψ_𝑖,0(𝑥,𝜇,𝜃) and Ψ_0,𝑗(𝑥,𝜇,𝜃) for 𝑖=1,…,𝐼 and 𝑗=1,…,𝐽 will satisfy one-dimensional transport problems subject to the same of boundary conditions of the original problem in the variable 𝑥.

3. Analysis

Now, we solve the first-order linear differential equation system with isotropic scattering, that is, 𝜎_𝑠(𝜇^′,𝜙^′→𝜇,𝜙)≡𝜎_𝑠= constant. Assuming isotropic scattering, the equation (2.13) is written as 𝜇𝜕Ψ_𝑖,𝑗𝜕𝑥𝑥,𝜇,𝜃+𝜎_𝑡Ψ_𝑖,𝑗𝑥,𝜇,𝜃=𝐺_𝑖,𝑗𝑥;𝜇,𝜃𝜎¹₋₁₀^2𝜋Ψ_𝑖,𝑗𝑥,𝜇^′,𝜃^′𝑑𝜃^′𝑑𝜇^′, (3.1) for 𝐱∈Ω∶={(𝑥,𝑦)∶0≤𝑥≤1,−1≤𝑦≤1}, 𝜇∈[−1,1], and 𝜃∈[0,2𝜋].

Then, we have the following theorem that is subject to the boundary conditions (2.14).

Theorem 3.1. Consider the integrodifferential equation (3.1) under the boundary conditions (2.14), then the function Ψ_𝑖,𝑗(𝑥,𝜇,𝜃) satisfy the following linear system of algebraic equations: ^𝑁_𝑛=0𝐷_𝑛,𝑚𝑝𝛽_{𝑛,𝑖,𝑗}𝑝,𝜃−𝜎_𝑠^𝑁_𝑛=0𝑝𝛼_{𝑛,𝑖,𝑗}𝑝,𝜃+𝜎_𝑡𝛼_{𝑛,𝑖,𝑗}=𝑝,𝜃₁₋₁𝐺_𝑖,𝑗𝑥,𝜇,𝜃𝑊^𝑒_𝑛𝜇𝑑𝜇+^𝑁_𝑛=0𝐷_𝑛,𝑚𝛽_{𝑛,𝑖,𝑗}(0,𝜃_𝑁𝑛=0),𝐷_𝑛,𝑚𝑝𝛼_{𝑛,𝑖,𝑗}𝑝,𝜃−𝜎_𝑠^𝑁_𝑛=0𝑝𝛽_𝑛,𝑘𝑝,𝜃+𝜎_𝑡𝛽_{𝑛,𝑖,𝑗}=𝑝,𝜃₁₋₁𝐺_𝑖,𝑗𝑥,𝜇,𝜃𝑊^𝑜_𝑛𝜇𝑑𝜇+^𝑁_𝑛=0𝐷_𝑛,𝑚𝛼_{𝑛,𝑖,𝑗}(0,𝜃). (3.2)

Proof. For this problem we expand the angular flux in terms of the Walsh function in the angular variable with its domain extended into the interval [−1,1]. To end this, the Walsh function 𝑊_𝑛(𝜇) are extended in an even and odd fashion as follows, see [14]: 𝑊^𝑒_𝑛𝜇=𝑊_𝑛𝜇𝑊if𝜇≥0,_𝑛𝑊−𝜇if𝜇<0,_𝑜𝑛𝜇=𝑊_𝑛𝜇if𝜇≥0,−𝑊_𝑛−𝜇if𝜇<0, (3.3) for 𝑛=0,1,…,𝑁. The important feature of this procedure relies on the fact that a function 𝑓(𝜇) is defined in the interval [−1,1] it can be expanded in terms of these extended functions in the manner: 𝑓𝜇=^∞_𝑛=0𝑎_𝑛𝑊^𝑒_𝑛𝜇+𝑏_𝑛𝑊^𝑜_𝑛𝜇, (3.4) where the coefficients 𝑎_𝑛 and 𝑏_𝑛 are determined as 𝑎_𝑛=12¹₋₁𝑓𝜇𝑊^𝑒_𝑛𝜇𝑏𝑑𝜇,_𝑛=12¹₋₁𝑓𝜇𝑊^𝑜_𝑛𝜇𝑑𝜇. (3.5) So, in order to use the Walsh function for the solution of the problem (3.1), the angular flux is approximated by the truncated expansion: Ψ_𝑖,𝑗𝑥,𝜇,𝜃=^𝑁_𝑛=0𝛼_{𝑛,𝑖,𝑗}(𝑥,𝜃)𝑊^𝑒_𝑛𝜇+𝛽_{𝑛,𝑖,𝑗}(𝑥,𝜃)𝑊^𝑜_𝑛𝜇. (3.6) Inserting this expansion into the linear transport (3.1), it turns out ^𝑁_𝑛=0𝜇𝜕𝛼_{𝑛,𝑖,𝑗}()𝜕𝑥𝑥,𝜃+𝜎_𝑡𝛼_{𝑛,𝑖,𝑗}(𝑥,𝜃)𝑊^𝑒_𝑛𝜇+𝜇𝜕𝛽_{𝑛,𝑖,𝑗}()𝜕𝑥𝑥,𝜃+𝜎_𝑡𝛽_{𝑛,𝑖,𝑗}(𝑥,𝜃)𝑊^𝑜_𝑛𝜇=^𝑁_𝑛=0𝜎_𝑠¹₋₁₀^2𝜋𝛼_{𝑛,𝑖,𝑗}𝑥,𝜃^′𝑊^𝑒_𝑛𝜇^′𝑑𝜃^′𝑑𝜇^′+¹₋₁₀^2𝜋𝛼_{𝑛,𝑖,𝑗}𝑥,𝜃^′𝑊^𝑜_𝑛𝜇^′𝑑𝜃^′𝑑𝜇^′+𝐺_𝑖,𝑗𝑥,𝜇,𝜃. (3.7) Multiplying (3.7) by 𝑊^𝑒_𝑚, 𝑚=0,…,𝑁 and integrating over the interval [−1,1], results: ^𝑁_𝑛=0𝜕𝛽_{𝑛,𝑖,𝑗}()𝜕𝑥𝑥,𝜃¹₋₁𝜇𝑊^𝑜_𝑛𝜇𝑊^𝑒_𝑛𝜇𝑑𝜇+𝜎_𝑡𝛼_{𝑛,𝑖,𝑗}(𝑥,𝜃)¹₋₁𝑊^𝑒_𝑛𝜇𝑊^𝑒_𝑚𝜇=𝑑𝜇_𝑁𝑛=0𝜎_𝑠0^2𝜋𝛼_{𝑛,𝑖,𝑗}𝑥,𝜃^′𝑑𝜃^′1₋₁𝑊^𝑜_𝑛𝜇^′𝑊^𝑜_𝑛𝜇^′𝑑𝜇^′+¹₋₁𝐺_𝑖,𝑗𝑥,𝜇,𝜃𝑊^𝑒_𝑛𝜇𝑑𝜇. (3.8) Similarly, multiplying (3.7) by 𝑊⁰_𝑚, 𝑚=0,…,𝑁 and integrating yields: ^𝑁_𝑛=0𝜕𝛼_{𝑛,𝑖,𝑗}()𝜕𝑥𝑥,𝜃¹₋₁𝜇𝑊^𝑜_𝑛𝜇𝑊^𝑒_𝑛𝜇𝑑𝜇+𝜎_𝑡𝛽_{𝑛,𝑖,𝑗}(𝑥,𝜃)¹₋₁𝑊⁰_𝑛𝜇𝑊⁰_𝑚𝜇=𝑑𝜇_𝑁𝑛=0𝜎_𝑠0^2𝜋𝛽_{𝑛,𝑖,𝑗}𝑥,𝜃^′𝑑𝜃^′1₋₁𝑊^𝑜_𝑛𝜇^′𝑊^𝑜_𝑛𝜇^′𝑑𝜇^′+¹₋₁𝐺_𝑖,𝑗𝑥,𝜇,𝜃𝑊⁰_𝑛𝜇𝑑𝜇. (3.9) The integrals appearing in (3.8) and (3.9) are known and are given [14] as 𝐷_𝑛,𝑚=12¹₋₁𝜇𝑊^𝑜_𝑛𝜇𝑊^𝑒_𝑚𝜇𝑑𝜇=¹₀𝜇𝑊_{(𝑛+𝑚)mod2}𝜇 (3.10) or 𝐷_𝑛,𝑚=12if𝑛=𝑚,−2^−(𝑘+2),if(𝑛+𝑚)mod2=2^𝑘𝑘natural,0atanothercase, (3.11) where the notation (𝑛+𝑚)mod2 denotes the mod2 sum of the binary digits 𝑛 and 𝑚 [15] ^𝑁_𝑛=0𝐷_𝑛,𝑚𝑝𝛽_{𝑛,𝑖,𝑗}𝑝,𝜃−𝜎_𝑠^𝑁_𝑛=0𝑝𝛼_{𝑛,𝑖,𝑗}𝑝,𝜃+𝜎_𝑡𝛼_{𝑛,𝑖,𝑗}=𝑝,𝜃₁₋₁𝐺_𝑖,𝑗𝑥,𝜇,𝜃𝑊^𝑒_𝑛𝜇𝑑𝜇+^𝑁_𝑛=0𝐷_𝑛,𝑚𝛽_{𝑛,𝑖,𝑗}(0,𝜃_𝑁𝑛=0),𝐷_𝑛,𝑚𝑝𝛼_{𝑛,𝑖,𝑗}𝑝,𝜃−𝜎_𝑠^𝑁_𝑛=0𝑝𝛽_𝑛,𝑘𝑝,𝜃+𝜎_𝑡𝛽_{𝑛,𝑖,𝑗}=𝑝,𝜃₁₋₁𝐺_𝑖,𝑗𝑥,𝜇,𝜃𝑊^𝑜_𝑛𝜇𝑑𝜇+^𝑁_𝑛=0𝐷_𝑛,𝑚𝛼_{𝑛,𝑖,𝑗}(0,𝜃). (3.12)

4. Operational Tau Method and Converting Transport Equation

The operational approach to the Tau method proposed by [16] describes converting of a given integral, integrodifferential equation or system of these equations to a system of linear algebraic equations based on three simple matrices: 0001𝛾=010001001⋯∶∶∶,𝜂=10020⋯0030⋯∶∶∶,𝑙=0100⋯2010⋯3⋯0⋯⋯∶∶∶. (4.1) We recall the following properties [17].

Lemma 4.1. Let 𝑢_𝑛(𝑥) be a polynomial as 𝑢_𝑛(𝑥)=^𝑛_𝑖=0𝑎_𝑖𝑥^𝑖=^∞_𝑖=0𝑎_𝑖𝑥^𝑖=𝑎_𝑛𝑋, (4.2) then we have(i)(𝑑^𝑟/𝑑𝑥^𝑟)𝑢_𝑛(𝑥)=𝑎_𝑛𝜂^𝑟𝑋, 𝑟=1,2,3,…,(ii)𝑥^𝑟𝑢_𝑛(𝑥)=𝑎_𝑛𝛾^𝑟𝑋, 𝑟=1,2,3,…,(iii)∫^𝑥_𝑎𝑢_𝑛(𝑥)𝑑𝑥=𝑎_𝑛𝑙𝑋|^𝑥_𝑎=𝑎_𝑛𝑙𝑋−𝑎_𝑛𝑙𝐴, where 𝑎_𝑛=[𝑎₀,𝑎₁,…,𝑎_𝑛,0,0,…] and 𝐴=[1,𝑎,𝑎²,…,𝑎^𝑛,…]^𝑇.

Now we solve (3.12) using the Tau method, for this we consider this equation because is equivalent 𝜇𝜕Ψ_𝑖,𝑗𝜕𝑥𝑥,𝜇,𝜃+𝜎_𝑡Ψ_𝑖,𝑗𝑥,𝜇,𝜃=𝐺_𝑖,𝑗𝑥;𝜇,𝜃+¹₋₁₀^2𝜋𝜎_𝑠𝜇^′,𝜃^′Ψ_𝑖,𝑗𝑥,𝜇^′,𝜃^′𝑑𝜃^′𝑑𝜇^′ (4.3) for 𝐱∈Ω∶={(𝑥,𝑦)∶0≤𝑥≤1,−1≤𝑦≤1}, 𝜇∈[−1,1], and 𝜃∈[0,2𝜋] subject to the following boundary conditions (2.14).

In order to convert (5.1) to a system of linear algebraic equations we define the linear differential operator 𝐷 of order 𝜁 with polynomial coefficients defined by 𝐷∶=^𝜁_𝑟=0𝑝_𝑟(𝑥)𝜕^𝑟𝜕𝑥^𝑟 (4.4) we will write for𝑝_𝑟(𝑥)∶=^𝛼𝑟_𝑘=0𝑝_𝑟𝑘𝑥^𝑘, (4.5) where 𝛼_𝑟 is the degree of 𝑝_𝑟(𝑥) and 𝑝_𝑟=(𝑝_𝑟0,…,𝑝_𝑟𝛼𝑟,0,0,…), Υ=(1,𝜇,𝜇²,…)^𝑇. We notice that 𝜇 is the independent variable and will be defined in a finite interval.

4.1. Matrix Representation for the Different Parts

Let 𝑉={𝑣₀(𝜇),𝑣₁(𝜇),…} be a polynomial basis given by 𝑉∶=𝑉Υ, where 𝑉 is nonsingular lower triangular matrix and degree (𝑣_𝑖(𝜇))≤𝑖, for 𝑖=0,1,2,…. Also for any matrix 𝑃, 𝑃_𝑣=𝑉𝑃𝑉⁻¹.

Matrix Representation for 𝐷Ψ_𝑖,𝑗(𝑥,𝜇,𝜃)
Theorem 4.2 (see [16]). For any linear differential operator 𝐷 defined by (5.2) and any series Ψ_𝑖,𝑗𝑥,𝜇,𝜃∶=𝑎𝑉,𝑎∶=Ψ_0,𝑗𝑥,𝜃^′,Ψ_1,𝑗𝑥,𝜃^′,Ψ_2,𝑗𝑥,𝜃^′,…,Ψ_𝑖,0𝑥,𝜃^′,Ψ_𝑖,1𝑥,𝜃^′,Ψ_𝑖,2𝑥,𝜃^′,…,, (4.6) we have 𝐷Ψ_𝑖,𝑗𝑥,𝜇,𝜃=𝑎Ξ_𝑣𝑉, (4.7) where Ξ_𝑣=𝑉^𝜁_𝑖=0𝜂^𝑖𝑝_𝑖𝛾𝑉⁻¹. (4.8)

Matrix Representation for the Integral Form
Let us assume that 𝜎_𝑠𝜇^′,𝜃^′=^𝑛_𝑗=0𝜎_𝑠𝜃^′𝑣_𝑗𝜇^′,Ψ_𝑖,𝑗𝑥,𝜇^′,𝜃^′=^𝑛_𝑖,𝑗=0Ψ_𝑖,𝑗𝑥,𝜃^′𝑣_𝑖𝜇^′=𝑎𝑉, (4.9) then we can write ¹₋₁₀^2𝜋𝜎_𝑠𝜇^′,𝜃^′Ψ_𝑖,𝑗𝑥,𝜇^′,𝜃^′𝑑𝜃^′𝑑𝜇^′=^𝑛_𝑖,𝑗=00^2𝜋𝜎_𝑠𝜃^′Ψ_𝑖,𝑗𝑥,𝜇^′,𝜃^′𝑑𝜃^′1₋₁𝑣_𝑗𝜇^′𝑣_𝑖𝜇^′𝑑𝜇^′=𝑎𝐾𝑉, (4.10) where 𝐾=^𝑛𝑗=0𝑘_0,𝑗𝜛_𝑗0⋯^𝑛_𝑖,𝑗=0𝑘_0,𝑗𝜛_𝑖𝑗00⋯⋯⋯⋯⋯⋯⋯_𝑛𝑗=0𝑘_0,𝑗𝜛_𝑖𝑗⋯^𝑛_𝑖,𝑗=0𝑘_0,𝑗𝜛_𝑖𝑗00⋯⋯⋯⋯⋯⋯⋯, (4.11) with 𝑘_𝑖,𝑗=∫₀^2𝜋𝜎_𝑠(𝜃^′), Ψ_𝑖,𝑗(𝑥,𝜇^′,𝜃^′)𝑑𝜃^′, and 𝜛_𝑖𝑗=∫¹₋₁𝑣_𝑖(𝜇^′)𝑣_𝑗(𝜇^′)𝑑𝜇^′ for 𝑖,𝑗=0,1,…,𝑛.

Matrix Representation for the Boundary Conditions
Ψ_𝑖,𝑗𝑥,𝜇^′,𝜃^′=^∞𝜁_𝑖=0𝑘=0Ψ_𝑗𝑘𝑥,𝜃^′𝑣_𝑗𝜇^′=𝑎ColΨ_0,𝑗𝑥,𝜃^′,Ψ_1,𝑗𝑥,𝜃^′,…,Ψ_𝑖,0𝑥,𝜃^′,Ψ_𝑖,1𝑥,𝜃^′,…,=𝑎𝐵_𝑗, (4.12) where 𝑗=1,…,𝜁. It follows from (5.4) and (5.5) that 𝜇𝜕Ψ_𝑖,𝑗𝜕𝑥𝑥,𝜇,𝜃+𝜎_𝑡Ψ_𝑖,𝑗𝑥,𝜇,𝜃=¹₋₁₀^2𝜋𝜎_𝑠𝜇^′,𝜃^′Ψ_𝑖,𝑗𝑥,𝜇^′,𝜃^′𝑑𝜃^′𝑑𝜇^′=1𝜇𝑎𝐾𝑉−𝑎Ξ_𝑣𝑉. (4.13) Let 𝐺_𝑖,𝑗(𝑥;𝜇,𝜃)=∑^𝑛_𝑑=0𝐺_{𝑖,𝑗𝑑}(𝑥;𝜃)𝑣_𝑑(𝜇)=𝐺_𝑖,𝑗(𝑥;𝜃)𝑉 with 𝐺_𝑖,𝑗=(𝐺_𝑖,𝑗0,…,𝐺_{𝑖,𝑗𝑛},0,0,…).

5. Error Estimation

Consider now the discrete ordinates (𝑆_𝑁) approximation of the equation (2.13) for 𝑚=1,…,𝑀,𝜇_𝑚𝜕Ψ_{𝛼,𝛽,𝑁}𝜕𝑥𝑥,𝜇_𝑚,𝜃+𝜎_𝑡Ψ_{𝛼,𝛽,𝑁}𝑥,𝜇_𝑚,𝜃=^𝑀_𝑛=1𝜔_𝑛Ψ_{𝛼,𝛽,𝑁}𝑥,𝜇_𝑚,𝜃+𝐺_𝛼,𝛽𝑥,𝜇_𝑚,𝜃. (5.1) In this part, we evaluate an error estimator for the approximate solution of (2.13), we suppose that equations (2.13) and (5.1) have the same boundary conditions. Let us call𝜖_𝑚𝑥,𝜇,𝜃=Ψ_𝑖,𝑗𝑥,𝜇_𝑚,𝜃−Ψ_{𝑖,𝑗,𝑚}𝑥,𝜇,𝜃 (5.2) this error function of the approximate solution Ψ_{𝑖,𝑗,𝑚} to Ψ_𝑖,𝑗 where Ψ_𝑖,𝑗 is the exact solution of (2.13). Hence, Ψ_𝑖,𝑗(𝑥,𝜇_𝑚,𝜃) satisfies the following problem:𝜇_𝑚𝜕Ψ_{𝛼,𝛽,𝑁}𝜕𝑥𝑥,𝜇_𝑚,𝜃+𝜎_𝑡Ψ_{𝛼,𝛽,𝑁}𝑥,𝜇_𝑚,𝜃=^𝑀_𝑛=1𝜔_𝑛Ψ_{𝛼,𝛽,𝑁}𝑥,𝜇_𝑚,𝜃+𝐺_𝛼,𝛽𝑥,𝜇_𝑚,𝜃+𝐻_𝑁𝜇_𝑚. (5.3) We can evaluate the perturbation term 𝐻_𝑁(𝜇_𝑚) by substituting the computed solution into the equation𝐻_𝑁𝜇_𝑚=𝜇_𝑚𝜕Ψ_{𝛼,𝛽,𝑁}𝜕𝑥𝑥,𝜇,𝜃+𝜎_𝑡Ψ_{𝛼,𝛽,𝑁}𝑥,𝜇,𝜃−^𝑀_𝑛=1𝜔_𝑛Ψ_{𝛼,𝛽,𝑁}𝑥,𝜇_𝑚,𝜃−𝐺_𝛼,𝛽𝑥,𝜇_𝑚,𝜃 (5.4) after doing some algebraic manipulations, the error functions 𝜖_𝑚(𝜇) satisfies the problem𝐻_𝑁𝜇_𝑚=𝜇_𝑚𝜕𝜖_𝑚𝑥,𝜇,𝜃𝜕𝑥+𝜎_𝑡𝜖_𝑚𝑥,𝜇,𝜃−^𝑀_𝑛=1𝜔_𝑛𝜖_𝑚𝑥,𝜇,𝜃. (5.5)

6. Conclusion

In general, obtaining solutions of some integrodifferential equations are usually difficult. In our recent works we have used Walsh functions, Chebyshev polynomials and Lengendre polynomials in order to reduces these kind of equations. However our present work suggests that the Tau method can be a good approximation to the exact solutions. The application of the Tau method by using the orthogonal polynomials will be considered as a future work.

Acknowledgments

The authors thank the referee(s) for the helpful and significant comments that bring the attention of authors to the references [7–11] which update the bibliography. The authors gratefully also acknowledge that this research was partially supported by the University Putra Malaysia under the ScienceFund Grant no. 06-01-04-SF1050 and Research University Grant Scheme no. 05-01-09-0720RU.