# VIM for Solving the Pollution Problem of a System of Lakes

Authors: J. Biazar; M. Shahbala; H. Ebrahimi
Journal: Journal of Control Science and Engineering
 (2010)

---

## Abstract

Pollution has become a very serious threat to our environment. Monitoring pollution is the first step toward planning to save the environment. The use of differential equations of monitoring pollution has become possible. In this paper the pollution problem of three lakes with interconnecting channels has been studied. The variational iteration method has been applied to compute an approximate solution of the system of differential equations, governing on the problem. Three different types of input models: sinusoidal, impulse, and step will be considered for monitoring the pollution in the lakes. The results are compared with those obtained by Adomian decomposition method. This comparison reveals that the variational iteration method is easier to be implemented.


---

## Body

## 1. Introduction

He’s variational iteration method plays an important role in recent researches in different fields of sciences and engineering. This method is proposed by the Chinese mathematician, He, [1–3] as a modification of the general Lagrange multiplier method. It has been shown to solve efficiently, easily, and accurately a large class of nonlinear problems with approximation converging to accurate solutions. This method has been used to study many problems in different disciplines, hyperbolic equation [4], wave equation [5], nonlinear chemistry problems [6], and quadratic Riccati differential equation [7].To illustrate its basic idea of the method, we consider the following nonlinear equation:(1)Lu(t)+Nu(t)=g(t),
where L is a linear operator, Nis a nonlinear operator, and g(t) is a known analytical function.According to VIM, we can construct the following correction functional:(2)un+1(t)=un(t)+∫0tλ{Lun(s)+Nũn(s)-g(s)}ds,
where λ is a Lagrange multiplier [8], which can be identified optimally via variational theory [9, 10], un is the nth approximate solution, and  ũn denotes a restricted variation [10], that is, δũn=0.


## 2. Modeling of the Pollution Problem

Figure1 shows the system of three lakes that are modeled in [11]. Each lake is considered to be a large compartment and the interconnecting channel as pipes between the compartments.Figure 1
System of three lakes with interconnecting channels.The direction of flow in the channels or pipes is indicated by the arrows in the figure. A pollutant is introduced into the first lake andp(t) denotes the rate at which the pollutant enters the lake per unit time. The function p(t) may be constant or may vary with time. We are interested in knowing the level of pollution in each lake at any time.Letxi(t) denotes the amount of the pollution in lake i at any time t≥0, where i=1,2,3. We assume the pollutant in each lake to be uniformly distributed throughout the lake by some mixing process, and the volume of water Vi in lake i remains constant for each of the lakes. Also we assume that the type of pollution is persistent and not degrading to other forms. Then the concentration of the pollutant in lake i at any time is given by
(3)ci(t)=xi(t)Vi.
Each lake initially is assumed to be free of any contaminant, so xi(0)=0 for each i=1,2,3. To model the dynamic behavior of the system of lakes, we let constant Fji denote the flow rate from lake i to lake j. The flux of pollutant flowing from lake i into lake j at any time t, denoted by rji(t), is defined by
(4)rji(t)=Fjici(t)=FjiVixi(t).
Thus rji(t) measures the rate at which the concentration of pollutant in lake i flows into lake j at time t. We will observe that
(5)rateofchangeofpollutant=Inputrate-outputrate.
Applying this principle to each lake results in the following system of first-order equations:
(6)dx1dt=F13V3x3(t)+p(t)-F31V1x1(t)-F21V1x1(t),dx2dt=F21V1x1(t)-F32V2x2(t),dx3dt=F31V1x1(t)+F32V2x2(t)-F13V3x3(t).
To remain constant the volume of each lake, the flow rate into each lake must balance the flow out of the lake. So we earn the following conditions:
(7)Lake1:F13=F21+F31,Lake2:F21=F32,Lake3:F31+F32=F13.


## 3. Numerical Results and Comparison with the ADM

There are three different types of input models used to monitor pollutants in a lake. The input models are the sinusoidal, impulse, and step input.


### 3.1. Sinusoidal Input

The Sinusoidal input model is used for pollutants that are introduced to the lake periodically. Sinusoidal input changesp(t) to a more useful one
(8)p(t)=pi*(1+a*sin(2π*tT)).The additional variables for the sinusoidal model are defined as follows:a:
the normalize amplitude,T:
the period of fluctuations,pi:
the average input concentration of pollutant.The actual amplitude of the term ispi(1+a) for the peaks and pi(1-a) for the valleys.An example would be a manufacturing plant dumping waste producing more output during the day than at night because of hour of operation, hence a periodic input. The concentration of the lake eventually converges to the average input concentration of the contaminant.Now we assumep(t)=1+sin(t), and for comparison with the results found by Biazar [11], the following values, for variables, are considered:
(9)V1=2900mi3,V2=850mi3,V3=1180mi3,F21=18mi3/year,F32=18mi3/year,F31=20mi3/year,F13=38mi3/year.So system (6) will be
(10)x1′(t)=381180x3(t)+(1+sin(t))-382900x1(t),x2′(t)=182900x1(t)-18850x2(t),x3′(t)=202900x1(t)+18850x2(t)-381180x3(t),x1(0)=0,x2(0)=0,x3(0)=0.According to VIM, we can construct the correction functionals of system (10) as follows:
(11)x1n+1(t)=x1n(t)+∫0tλ1{x1n′(ξ)-381180x̃3n(ξ)-1-sin(ξ)+382900x̃1n(ξ)}dξ,x2n+1(t)=x2n(t)+∫0tλ2{x2n′(ξ)-182900x̃1n(ξ)+18850x̃2n(ξ)}dξ,x3n+1(t)=x3n(t)+∫0tλ3{x3n′(ξ)-202900x̃1n(ξ)-18850x̃2n(ξ)+381180x̃3n(ξ)}dξ,
where λ1, λ2, and λ3 are general Lagrange multipliers, x̃1n(ξ), x̃2n(ξ), and x̃3n(ξ) denote restricted variations, that is, δx̃1n(ξ)=δx̃2n(ξ)=δx̃3n(ξ)=0.Making the above correction functional stationary, we have(12)δx1n+1(t)=δx1n(t)+λ1(ξ)δx1n(ξ)|0t-∫0tλ1′(ξ)δx1n(ξ)dξ=0,δx2n+1(t)=δx2n(t)+λ2(ξ)δx2n(ξ)|0t-∫0tλ2′(ξ)δx2n(ξ)dξ=0,δx3n+1(t)=δx3n(t)+λ3(ξ)δx3n(ξ)|0t-∫0tλ3′(ξ)δx3n(ξ)dξ=0.The stationary conditions can be obtained as follows:(13)δx1n:{1+λ1(ξ)|ξ=t=0,λ1′(ξ)=0,δx2n:{1+λ2(ξ)|ξ=t=0,λ2′(ξ)=0,δx3n:{1+λ3(ξ)|ξ=t=0,λ3′(ξ)=0.Therefore, the Lagrange multipliers can be identified as(14)λ1=λ2=λ3=-1.Substituting (14) into the correction functionals (11) results in the following iteration formula:
(15)x1n+1(t)=x1n(t)-∫0t{x1n′(ξ)-381180x3n(ξ)-(1+sin(ξ))+382900x1n(ξ)}dξ,x2n+1(t)=x2n(t)-∫0t{x2n′(ξ)-182900x1n(ξ)+18850x2n(ξ)}dξ,x3n+1(t)=x3n(t)-∫0t{x3n′(ξ)-202900x1n(ξ)-18850x2n(ξ)+381180x3n(ξ)}dξ.We start with initial approximationsx10=x1(0)=0,x20=x2(0)=0,x30=x3(0)=0.By the above iteration formula, we can obtain a few first terms being calculated:(16)x11=1+t-cos(t),x21=0,x31=0,x12=1+14311450t-cos(t)-192900t2+191450sin(t),x22=91450t+92900t2-91450sin(t),x32=1145t+1290t2-1145sin(t),x13=123998651124047500+14311450t-123998651124047500cos(t)-1576601248095000t2+191450sin(t)+16283248095000t3,x23=91450t+4284914297000t2-91450sin(t)+15217148500-50714297000t3-15217148500cos(t),x33=1145t+70808812108807500t2-1145sin(t)+1908691054403750-636232108807500t3-1908691054403750cos(t).⋮Continuing in this manner, we can find the rest of components.A fewm terms approximation to the solutions are considered as
(17)x1(t)≈x1m,x2(t)≈x2m,x3(t)≈x3m.Some numerical values of solutions are presented in Table1 with n being the order of iteration formula of VIM.Table 1
Result for sinusoidal model.
ntx1(t)x2(t)x3(t)10.000.00000000.00000000.000000020.010.1004934× 10-10.3113793× 10-60.3459770× 10-630.020.2019735× 10-10.1249370× 10-50.1388263× 10-540.030.3044401× 10-10.2820067× 10-50.3133661× 10-550.040.4078927× 10-10.5029427× 10-50.5588849× 10-560.050.5123310× 10-10.7883429× 10-50.8760533× 10-570.060.6177546× 10-10.1138805× 10-40.1265541× 10-480.070.7241617× 10-10.1554927× 10-40.1728018× 10-490.080.8315528× 10-10.2037305× 10-40.2264153× 10-4100.090.9399265× 10-10.2586535× 10-40.2874615× 10-4110.100.10492820.3203213× 10-40.3560070× 10-4


### 3.2. Impulse Input

The impulse input model is used for pollutants that have been released into the lake immediately. Impulse input functions have a spike, and everywhere else the function is zero. The spike represents the time at which the contaminant was dumped.An example would be dumping one barrel of gasoline into a lake at time zero. In one short heavy dumping of waste, with no reaction with the lake, inputp(t)  is equal to 100 units for a duration of 10 units of time.Now we assumep(t)=100, and so the system (6) with parameters that are defined in Section 3.1 will be
(18)x1′(t)=381180x3(t)+100-382900x1(t),x2′(t)=182900x1(t)-18850x2(t),x3′(t)=202900x1(t)+18850x2(t)-381180x3(t),x1(0)=0,x2(0)=0,x3(0)=0.Therefore, by solving system (18) with VIM, as have been shown in Sinusoidal input model, the following results are earned:
(19)x11=100t,x21=0,x31=0,x12=100t-1929t2,x22=929t2,x32=1029t2,x13=100t-1929t2+162832480950t3,x23=929t2-507142970t3,x33=1029t2-6362321088075t3.⋮Some numerical values of solutions are presented in Table2 with n being the order of iteration formula of VIM.Table 2
Result for impulse model.
ntx1(t)x2(t)x3(t)10.000.00000000.00000000.000000020.010.99893450.3103448× 10-40.3448275× 10-430.020.1999738× 1010.1241095× 10-30.1379069× 10-340.030.2999411× 1010.2792146× 10-30.3102634× 10-350.040.3998952× 1010.4963248× 10-30.5515311× 10-360.050.4998363× 1010.7754190× 10-30.8616919× 10-370.060.5997643× 1010.1116476× 10-20.1240728× 10-280.070.6996792× 1010.1519474× 10-20.1688621× 10-290.080.7995810× 1010.1984392× 10-20.2205353× 10-2100.090.8994698× 1010.2511210× 10-20.2790905× 10-2110.100.9993455× 1010.3099905× 10-20.3445261× 10-2


### 3.3. Step Model

The step input model is used for pollutants that enter the lake at steady concentration and rate and continue the same way indefinitely. The pollutant enters the system at time zero, and the concentration previous to time zero was nothing. The significant points for step input are that the input increases suddenly at time zero and that the input remains essentially constant after time zero.An example would be when a manufacturing plant starts production at time zero and starts dumping raw sewage at a constant rate and concentration.Now we assumep(t)=100t, so system (6) with parameters that are defined in Section 3.1 would be
(20)x1′(t)=381180x3(t)+100t-382900x1(t),x2′(t)=182900x1(t)-18850x2(t),x3′(t)=202900x1(t)+18850x2(t)-381180x3(t),x1(0)=0,x2(0)=0,x3(0)=0.Now by using VIM for system (20), the following approximations for the solution are earned:
(21)x11=50t2,x21=0,x31=0,x12=50t2-1987t3,x22=329t3,x32=1087t3,x13=50t2-1987t3+162839923800t4,x23=329t3-507571880t4,x33=1087t3-6362384352300t4.⋮Some numerical values of solutions are presented in Table3 with n being the order of iteration formula of VIM.Table 3
Result for step model.
ntx1(t)x2(t)x3(t)10.000.00000000.00000000.000000020.010.4999782× 10-20.1034448× 10-60.1149425× 10-630.020.1999825× 10-10.8274444× 10-50.9194195× 10-640.030.4499410× 10-10.2792385× 10-50.3102837× 10-550.040.7998603× 10-10.6618421× 10-50.7354391× 10-560.050.12497270.1292549× 10-40.1436103× 10-470.060.17995280.2233334× 10-40.2481781× 10-480.070.24492510.3546148× 10-40.3940718× 10-490.080.31988820.5292922× 10-40.5881969× 10-4100.090.40484090.7535566× 10-40.8374364× 10-4110.100.49978180.1033597× 10-30.1148671× 10-3


## 3.1. Sinusoidal Input

The Sinusoidal input model is used for pollutants that are introduced to the lake periodically. Sinusoidal input changesp(t) to a more useful one
(8)p(t)=pi*(1+a*sin(2π*tT)).The additional variables for the sinusoidal model are defined as follows:a:
the normalize amplitude,T:
the period of fluctuations,pi:
the average input concentration of pollutant.The actual amplitude of the term ispi(1+a) for the peaks and pi(1-a) for the valleys.An example would be a manufacturing plant dumping waste producing more output during the day than at night because of hour of operation, hence a periodic input. The concentration of the lake eventually converges to the average input concentration of the contaminant.Now we assumep(t)=1+sin(t), and for comparison with the results found by Biazar [11], the following values, for variables, are considered:
(9)V1=2900mi3,V2=850mi3,V3=1180mi3,F21=18mi3/year,F32=18mi3/year,F31=20mi3/year,F13=38mi3/year.So system (6) will be
(10)x1′(t)=381180x3(t)+(1+sin(t))-382900x1(t),x2′(t)=182900x1(t)-18850x2(t),x3′(t)=202900x1(t)+18850x2(t)-381180x3(t),x1(0)=0,x2(0)=0,x3(0)=0.According to VIM, we can construct the correction functionals of system (10) as follows:
(11)x1n+1(t)=x1n(t)+∫0tλ1{x1n′(ξ)-381180x̃3n(ξ)-1-sin(ξ)+382900x̃1n(ξ)}dξ,x2n+1(t)=x2n(t)+∫0tλ2{x2n′(ξ)-182900x̃1n(ξ)+18850x̃2n(ξ)}dξ,x3n+1(t)=x3n(t)+∫0tλ3{x3n′(ξ)-202900x̃1n(ξ)-18850x̃2n(ξ)+381180x̃3n(ξ)}dξ,
where λ1, λ2, and λ3 are general Lagrange multipliers, x̃1n(ξ), x̃2n(ξ), and x̃3n(ξ) denote restricted variations, that is, δx̃1n(ξ)=δx̃2n(ξ)=δx̃3n(ξ)=0.Making the above correction functional stationary, we have(12)δx1n+1(t)=δx1n(t)+λ1(ξ)δx1n(ξ)|0t-∫0tλ1′(ξ)δx1n(ξ)dξ=0,δx2n+1(t)=δx2n(t)+λ2(ξ)δx2n(ξ)|0t-∫0tλ2′(ξ)δx2n(ξ)dξ=0,δx3n+1(t)=δx3n(t)+λ3(ξ)δx3n(ξ)|0t-∫0tλ3′(ξ)δx3n(ξ)dξ=0.The stationary conditions can be obtained as follows:(13)δx1n:{1+λ1(ξ)|ξ=t=0,λ1′(ξ)=0,δx2n:{1+λ2(ξ)|ξ=t=0,λ2′(ξ)=0,δx3n:{1+λ3(ξ)|ξ=t=0,λ3′(ξ)=0.Therefore, the Lagrange multipliers can be identified as(14)λ1=λ2=λ3=-1.Substituting (14) into the correction functionals (11) results in the following iteration formula:
(15)x1n+1(t)=x1n(t)-∫0t{x1n′(ξ)-381180x3n(ξ)-(1+sin(ξ))+382900x1n(ξ)}dξ,x2n+1(t)=x2n(t)-∫0t{x2n′(ξ)-182900x1n(ξ)+18850x2n(ξ)}dξ,x3n+1(t)=x3n(t)-∫0t{x3n′(ξ)-202900x1n(ξ)-18850x2n(ξ)+381180x3n(ξ)}dξ.We start with initial approximationsx10=x1(0)=0,x20=x2(0)=0,x30=x3(0)=0.By the above iteration formula, we can obtain a few first terms being calculated:(16)x11=1+t-cos(t),x21=0,x31=0,x12=1+14311450t-cos(t)-192900t2+191450sin(t),x22=91450t+92900t2-91450sin(t),x32=1145t+1290t2-1145sin(t),x13=123998651124047500+14311450t-123998651124047500cos(t)-1576601248095000t2+191450sin(t)+16283248095000t3,x23=91450t+4284914297000t2-91450sin(t)+15217148500-50714297000t3-15217148500cos(t),x33=1145t+70808812108807500t2-1145sin(t)+1908691054403750-636232108807500t3-1908691054403750cos(t).⋮Continuing in this manner, we can find the rest of components.A fewm terms approximation to the solutions are considered as
(17)x1(t)≈x1m,x2(t)≈x2m,x3(t)≈x3m.Some numerical values of solutions are presented in Table1 with n being the order of iteration formula of VIM.Table 1
Result for sinusoidal model.
ntx1(t)x2(t)x3(t)10.000.00000000.00000000.000000020.010.1004934× 10-10.3113793× 10-60.3459770× 10-630.020.2019735× 10-10.1249370× 10-50.1388263× 10-540.030.3044401× 10-10.2820067× 10-50.3133661× 10-550.040.4078927× 10-10.5029427× 10-50.5588849× 10-560.050.5123310× 10-10.7883429× 10-50.8760533× 10-570.060.6177546× 10-10.1138805× 10-40.1265541× 10-480.070.7241617× 10-10.1554927× 10-40.1728018× 10-490.080.8315528× 10-10.2037305× 10-40.2264153× 10-4100.090.9399265× 10-10.2586535× 10-40.2874615× 10-4110.100.10492820.3203213× 10-40.3560070× 10-4


## 3.2. Impulse Input

The impulse input model is used for pollutants that have been released into the lake immediately. Impulse input functions have a spike, and everywhere else the function is zero. The spike represents the time at which the contaminant was dumped.An example would be dumping one barrel of gasoline into a lake at time zero. In one short heavy dumping of waste, with no reaction with the lake, inputp(t)  is equal to 100 units for a duration of 10 units of time.Now we assumep(t)=100, and so the system (6) with parameters that are defined in Section 3.1 will be
(18)x1′(t)=381180x3(t)+100-382900x1(t),x2′(t)=182900x1(t)-18850x2(t),x3′(t)=202900x1(t)+18850x2(t)-381180x3(t),x1(0)=0,x2(0)=0,x3(0)=0.Therefore, by solving system (18) with VIM, as have been shown in Sinusoidal input model, the following results are earned:
(19)x11=100t,x21=0,x31=0,x12=100t-1929t2,x22=929t2,x32=1029t2,x13=100t-1929t2+162832480950t3,x23=929t2-507142970t3,x33=1029t2-6362321088075t3.⋮Some numerical values of solutions are presented in Table2 with n being the order of iteration formula of VIM.Table 2
Result for impulse model.
ntx1(t)x2(t)x3(t)10.000.00000000.00000000.000000020.010.99893450.3103448× 10-40.3448275× 10-430.020.1999738× 1010.1241095× 10-30.1379069× 10-340.030.2999411× 1010.2792146× 10-30.3102634× 10-350.040.3998952× 1010.4963248× 10-30.5515311× 10-360.050.4998363× 1010.7754190× 10-30.8616919× 10-370.060.5997643× 1010.1116476× 10-20.1240728× 10-280.070.6996792× 1010.1519474× 10-20.1688621× 10-290.080.7995810× 1010.1984392× 10-20.2205353× 10-2100.090.8994698× 1010.2511210× 10-20.2790905× 10-2110.100.9993455× 1010.3099905× 10-20.3445261× 10-2


## 3.3. Step Model

The step input model is used for pollutants that enter the lake at steady concentration and rate and continue the same way indefinitely. The pollutant enters the system at time zero, and the concentration previous to time zero was nothing. The significant points for step input are that the input increases suddenly at time zero and that the input remains essentially constant after time zero.An example would be when a manufacturing plant starts production at time zero and starts dumping raw sewage at a constant rate and concentration.Now we assumep(t)=100t, so system (6) with parameters that are defined in Section 3.1 would be
(20)x1′(t)=381180x3(t)+100t-382900x1(t),x2′(t)=182900x1(t)-18850x2(t),x3′(t)=202900x1(t)+18850x2(t)-381180x3(t),x1(0)=0,x2(0)=0,x3(0)=0.Now by using VIM for system (20), the following approximations for the solution are earned:
(21)x11=50t2,x21=0,x31=0,x12=50t2-1987t3,x22=329t3,x32=1087t3,x13=50t2-1987t3+162839923800t4,x23=329t3-507571880t4,x33=1087t3-6362384352300t4.⋮Some numerical values of solutions are presented in Table3 with n being the order of iteration formula of VIM.Table 3
Result for step model.
ntx1(t)x2(t)x3(t)10.000.00000000.00000000.000000020.010.4999782× 10-20.1034448× 10-60.1149425× 10-630.020.1999825× 10-10.8274444× 10-50.9194195× 10-640.030.4499410× 10-10.2792385× 10-50.3102837× 10-550.040.7998603× 10-10.6618421× 10-50.7354391× 10-560.050.12497270.1292549× 10-40.1436103× 10-470.060.17995280.2233334× 10-40.2481781× 10-480.070.24492510.3546148× 10-40.3940718× 10-490.080.31988820.5292922× 10-40.5881969× 10-4100.090.40484090.7535566× 10-40.8374364× 10-4110.100.49978180.1033597× 10-30.1148671× 10-3


## 4. Conclusion and Discussion

In this paper, the problem of pollutions of three lakes with interconnecting channels has been considered. Different input models have been used for monitoring the pollution in three lakes. The variational iteration method as a powerful device has been used to solve the system of differential equations, governing on the problem. For different input models, the results of applyingn times iteration formula in VIM procedure for different times t (year) are shown in suitable tables. The results are compared with those obtained by ADM in [11]. This comparison shows that the results are the same, but the solution procedure of the VIM is simpler than of ADM and the amounts of computations required in VIM is much lesser than ADM.



---

*Source: 829152-2010-09-06.xml*