Solving Load Flow Problems of Power System by Explicit Pseudo-Transient Continuation ( E-Ψtc ) Method

This paper is a further study of two papers [1] and [2], which were related to Ill-Conditioned Load Flow Problems and were published by IEEE Trans. PAS. The authors of this paper have some different opinions, for example, the 11-bus system is not an ill-conditioned system. In addition, a new approach to solve Load Flow Problems, E-Ψtc, is introduced. It is an explicit method; solving linear equations is not needed. It can handle very tough and very large systems. The advantage of this method has been fully proved by two examples. The authors give this new method a detailed description of how to use it to solve Load Flow Problems and successfully apply it to the 43-bus and the 11-bus systems. The authors also propose a strategy to test the reliability, and by solving gradient equations, this new method can answer if the solution exists or not.


Introduction
In the early 1980's, there were two papers which studied ill-conditioned Load Flow Problems ( [1] and [2]).Different numerical method was used to analyze the ill-conditioned equations and got some useful results.Because of the importance of this problem in power system, it is necessary to study this problem again by using new numerical method.Since 1960's, people have taken the direct method as a major tool to solve load flow equations; among them, Newton-Gaussian elimination (LU) is a typical representation [3].Despite the fact that some people proposed nonlinear programming method [4], this method was not widely used.
Recently people have turned their attention to iterative method, e.g., Inexact-Newton-GMRES method, ([5] and [6]), but none of them is involved with the ill-conditioned problems.
The extent of ill-conditioned problem is determined by the conditioned number of the system.From [2], for 11-bus and 43-bus systems, they are 0.1086 × 10 4 and 0.2560 × 10 5 respectively.For double precision calculation, it is not a severe problem; we can determine the existence of the solutions of the two systems by the normal Newton method without any trouble.
For the conclusion (3) of [1] and the statements "Supposing that the correction vector Δx is obtained by some way…".For those statements, the authors of this paper have a different opinion.The details of discussion can be found in Appendix.
Brown method was used by [2].Ref. [7] did a deep and detailed study for this method and pointed out: "The method did not have a clear algorithmic description and hence was hard to implement".Maybe even more important is that for each iteration, it needs ( ) O n arithmetic operations.[8] improved Brown method, and also gave a clearer algorithmic description and reduced arithmetic operations from ( ) O n to ( ) O n .[2] did not take advantage of this point, and still used original Brown method.
Considering that the two papers were published in 1980's, we should not make excessive demands for authors of [2], but the conclusion for 11-bus and 43-bus systems, Newton-Raphson method is "divergent" (Table 2 on page 3653 of Ref. [2]), which is incorrect and conflicts with [1].

Description of the New Algorithm
Before the description of the new method, let us first review the derivation process of the two classical methods: 1) Newton method: From ( ) ( ) ( ) Or, for ( ) ( ) using explicit Euler numerical integration method: ( ) ( ) ( ) Take 1 h = , hence we have: That is to say, we can use both Static State and Dynamic State approaches to derive Newton method.Despite the Newton method can be recognized as an explicit ODE Euler method, the right side function includes an inverse matrix of ( ) It still belongs to implicit method.If the Newton method fails to converge, two approaches can be taken to improve the convergence: a) Using damp factor: Replacing x ∆ by ( )( ) < .This is equivalent to reducing the step size of Eu- ler method, making ( ) to fall in the stability region of Euler method.Here ρ is the spectral radius of the Jacobian matrix of ( ) ( ) b) Altering initial values: If the initial value ( ) closes to the solution x * , the Jacobian matrix of ( ) ( ) will close to unit matrix.Such that taking large step size h = 1, the stability of numerical integration can be kept.
2) About the fixed iteration For + , here h is the parameter.Fixed iteration method: ( ) will have the same solution of ( ) 0 h F x * = and ( ) 0 We can still interpret it as an explicit ODE Euler method.In fact, for differential equation ( ) Take h = 1, we can get classical fixed iteration: ( ) ( ) It is well known that the convergence condition of fixed iteration is: If h = 1, the convergence condition of classical Fixed Point Iteration Method does not satisfy.We can reduce h, if only all of the eigenvalues of the Jacobian matrix of ( ) We know that explicit ODE method only has very small stability region.If the Jacobian matrix of ( ) F x has large condition number, or the system is stiff, the conventional explicit ODE method will be very inefficient.
3) Does the explicit ODE method with the large stability region exist?The answer is positive.In [9] the authors of this paper proposed a new ODE method which has a very large stability region.Choosing proper parameter, the region in real axis direction can reach ∞ point.It seems strange, the textbook tells us the stability region of the explicit ODE method must be bounded, however our method includes a parameter ε.It plays a very important role for stability and [11] gives it a reasonable explanation.
Our purpose is to reach the steady point of the dynamic system ( ) as quickly as possible.We do not require the trajectory produced by our method approximates the solution curve of differential equation of ( ) Or we can say that our method is allowed to have large errors, if only the trajectory produced by our method is within the attracting domain of the steady point.So our method can take the very large step size, and goes to steady point quickly.The traditional explicit ODE methods cannot take large step size.Usually, their trajectory is very close to the solution of the system.In a sense the trajectory approximately reflects the real transient process of the system, however the trajectory produced by our method does not.So, our method can be recognized as "Pseudo-Transient Continuation." The existing "Pseudo-Transient Continuation" method can be written as follows: Various algorithms adopt different strategies to control h.They are known as the "LINEARLY IMPLICIT EULER DISCRETIZATION", [12].
Because each iteration step needs to solve a system of linear equations, so all of them are implicit methods.
Our method is the unique Explicit "Pseudo-Transient Continuation" method.In order to express this difference, our method is named by [11] as "Explicit Pseudo-Transient Continue" method.
From [9] and [13], we can write our scheme as follows: For ( ) 3) Here ( ) x n is the approximation of x * , the solution we expect.

The Advantage of the New Method
1) Our method is an explicit stiff ODE method, so it has a very strong capability to deal with tough problems.Brown's almost linear equation is as follows: [7] [14] ( ) T 0 1 0, for 1, 2, 1 1 0 0.5, 0.5, 0.5  is a solution.For n = 10, 30, Newton method diverges drastically, but Brown method converges well.If n > 35, according to [7], Brown method will not work.This has been proved by [14], for n = 40.What happens to our method?For n = 40, 100, 200, 400, we can get the solution easily.In this paper, we just show the result of n = 400: Virtually, the differential equation to deal with is: ( ) ( ) Here the diagonal matrix ( ) , the division has been avoided, or we put ( ) 1.0 n d x = .As we did in [9], we separated the calculation into three stages (in every stage we expect the 2-Norm of ( ) F x is descending about 5 magnitude order) and took parameters 2.0 400.0 ε = , Tol1 = 1D−0, Tol2 = 1D−5, Tol3 = 1D−10; h1 = 0.015, h2 = 0.15, h3 = 0.8.After 826 iterations, we got the solution: 1.00000000 Incidentally, the efficiency of Brown method is strongly related to the order of the equations.It is more efficient if the equations are ordered by increasing the degree of nonlinearity.In our method, each component can be calculated independently, so there is no ordering problem.
An experiment has been done: We place the nonlinear equation on the top of the ( ) F x , i.e. the first equation.Both of them have exactly the same results.
2) Our method is as simple as fixed iteration, so it can handle very large systems.Here we give an ill-conditioned, non-symmetric linear equations with dimension n = 1,000,000: .
We ask the solution is .
If E-Ψtc method is used, taking ε = 0.5.For any step size the numerical integration is stable.The only possible comparing method is fixed iteration (Euler) method.The maximum step size is h = 0.8.In order to get a good stability, we take h = 0.75.The initial Euclid norm is 0.7217D+03 After 150,000 fixed iterations, we have the following results: 3) Another outstanding feature is that our method is very suitable for parallel calculation.

Load Flow Calculation
Now we turn our focus to the main topic of this paper-Load Flow Calculation.
1) The Construction of the Differential Equations.The load flow equations can be formulated by a system of nonlinear.Equations: Here ( ) . In rectangular coordinate system, for PQ buses, the power equations are: For PV buses, the voltage magnitude values V is are specified and have the relation: and (3.3) constitute the power equations for PV buses.
For elimination method, the optimal ordering of the buses is the first task we need to solve.In our method, this problem does not need to be considered, because the memory of space only depends on the number of buses and not related to the ordering.For the calculating process, as we have pointed out in Section 3, in our method, each component can be calculated independently.
But for every bus, which corresponds to two real number equations.If the nonlinear equations , , , , , , , The corresponding differential equation should be arranged as follows: Because in the i th bus, the active power P i is mainly concerned with f i and reactive power Q i is mainly concerned with e i ([15] p. 205).The diagonal element of the Jacobian of ( ) F x will be the one which has the larg- est absolute value in each row.This is crucial for any iterative method.
In fact, for elimination method, every bus corresponds to a 2-dimension sub-block and which is also arranged just the same.
In order to let the readers understand the discussion above more clearly, let us give a further statement by a simple example: Consider a 2-dimensional linear equations: The matrix-vector form is: Here ( ) ( ) [ ] This expression is not correct.It will make the Jacobian matrix loose diagonal dominant.That is to say, for each row of the Jacobian matrix, the element with the largest absolute value will not be diagonal element.Obviously, for PV-bus, the one of diagonal elements will be "−2f".If flat start is used, then for the first step, the diagonal element will be zero.It will be the worst thing for any algorithm to solve the linear algebraic equations.For elimination method, almost all code, column pivot technique was taken.It automatically adjusts the order of equations so people usually ignore the "order problem".But for any iterative method, the "order problem" is crucial.
The correct expression should be: , , col. , x f e f e = .
(Ex.) ( ) ( ) ( ) ( ) ( ) Or the matrix form: For the details, please see [15] page 184.We can also change the order of equations; (Ex.) [ ] 2) Preconditioning We will give two preconditioning strategies: a) point diagonal preconditioning: The preconditioner is a diagonal matrix of the Jacobian.If you do not want to calculate the diagonal elements, according to the observation of the calculated data, we can further simplify this strategy.
For PQ bus: ( ) The precondition matrix can be simplified as a constant diagonal matrix.The elements are just B ii , and −2.0.It was named as D 1 .The differential equation that we deal with is: b) 2-dimension block diagonal preconditioning: We know every bus corresponds to a 2-dimension diagonal sub-block of the Jacobian matrix.These sub-blocks constitute a block diagonal matrix D 2 : The differential equation that we deal with is: ( ) ( ) The preconditioning not only raises efficiency, but also helps us to determine the parameter ε = 0.5.Because the Jacobian matrix diagonal elements of ( ) ( )

D F x
− are very close to 1, and it is fully possible, ∞-Norm of the Jacobian matrix < 3.0 (because the largest absolute value element of the Jacobian matrix is in the diagonal place of each row and the sparse property of the matrix) and ( ) 4 0.5 3 ρ * < , ρ is the spectral radius.(See the reference [10]).By this means, there is only one parameter h which needs be determined.
For ordinary load flow problem, if point diagonal preconditioning is used, h is in the range of 5 to 7 is suggested.If an iteration cannot converge, we suggest to reduce h.If h is small enough, say h < 1, the iteration still cannot converge, the iteration should be stopped, then check the whole process to see what mistake was taken.

43-Bus System
For this system, 8 calculating results are given, we take the Largest Absolute Values of the Mismatch (LAVM), less than 1D−05 as the stopping condition (double precision): 1) Normal Newton method: 7 Newton iteration, The LAVM = 0.2123D−08.
3) Euler method: step size h = 0.00065 (maximum step-length to meet the stability requirements), 157006 function evaluations.The LAVM = 1D−05.What conclusion that we can get based on the 1)-7) results?
From 3) we can see that this system is surely an ill-conditioned system, but from 1), the ill-conditioned extent is not severe, using double precision, only 7 Newton iterations are needed to make the LAVM = 0.2123D−08.
Compare 1) and 2), Optimal Multiplier Newton method saved one iteration, if the stopping condition is 1D−08.This is the contribution of [1].
Compare 3) and 7): the proportion of numbers of function evaluation is 157006/1335 ≈ 118.For 7), besides function evaluations, forty-two 2-dimension sub-matrix need to be inverted in each iteration.The workload of 2-dimension matrix inverse is relatively very small, so we can say that for ill-conditioned Load Flow Problems, using 7) is a good choice, especially when the system is large.
For 7), if someone considers that the workload of 1335 function evaluations is too big to be accepted, here we give another choice.This choice is based on the following observation: For Newton method, the quadratic convergence superiority can be utilized only when the norm of ( ) F x is small, say it is less than 0.1, but for func- tion iteration methods, at the beginning of the iteration, the norm has faster descending speed.So a hybrid calculation strategy is recommended.
Here, what we want to say is as following: For any existing Newton program, you can use our method as a tool to choose a good starting point.The textbook tells us: "The HOMOTOPY method can be used to generate a good starting value".So our method can be considered as a homotopy method.It can be used independently and can also be embedded into the Newton program as a tool to choose starting value.
There are a lot of references related to "HOMOTOPY", we merely list three of them [17]- [19].
For ordinary load flow problem, e.g.IEEE-39 IEEE-118, our method is also suitable.You can use the method by two ways: One is directly solve nonlinear equations, like what we did in this paper.[20], and another one is to solve the linear equations arising from the Newton method [21].
For nonlinear equations, the solution is not unique, as a research paper, we should give the final numerical results or some description for large systems, but neither did [1] nor [2].In this paper, listing all the 8 final results is not practical, because of the limited space.Here we only point out that no matter the LAVM is 1D-05 or 1D-12, the final results are almost the same.For bus 43, all of them are exactly the same: voltage magnitude is 1.0763, angle is −11.176.

11-Bus System
At the beginning, the initial LAVM = 0.1232E+01, we give the following results: 1) As [1] did, the normal Newton method, after 18 iterations, we have the results (Single precision): 0.3046E+00, 0. If we use double precision, we have almost the same results as above 1) and 2).Based on this observation, we can see: a) Oscillation occurred after 7 normal Newton method iterations, but for the optimal multiplier method, oscillation did not appear.
b) The LAVM cannot go down to zero and this is not caused by the defect of numerical precision of the computer and/or programming code.c) After 7 iterations, the LAVM of the optimal multiplier method is 0.2464E−03.After 7 iterations the LAVM of the normal Newton method is 0.2557E−03.
Unless it has a very small oscillation, the LAVM basically keeps monotone descending.(Both single and double precision have the same results).The LAVM is about 1/2 of the optimal multiplier method.4) When using Inexact Newton-GMRES (15), we use single and double precision to calculate, the results are almost the same.Here we give the double precision result: (For every GMRES loop, we calculate the initial residue R 0 and the final residue R 1 of the linear equations, and take R l /R 0 < 0.1, or Maximum Number of Restart = 100 as the stop condition for each Newton iteration.) That is to say, after 7 Newton iterations, there are totally (1 + 3 + 4 + 100 + 72 + 100 + 100 = 380) 380 × 15 = 5700 GMRES iterations, finally the LAVM reaches 0.2796D−02, then it no longer descends.
All results we got show that for the 11-bus system, the LAVM cannot go down to zero, in another word, for this system, the solution of ( ) 0 F x = , does not exist.What we should do is to seek the minimum of the func- tion

Solving Gradient Equation to Seek the Minimum for Bus-11
The first choice should be the Gauss-Newton method: if the Jacobi is a N N * nonsingular matrix, Gauss-Newton method will change into Newton method.As we have pointed out after 7 iterations the oscillation occurred.
We can also solve the gradient equations by using Newton method.In order to compare the two results, re-solve it again by LINPACK SUBROUTINE DGECO.By this way, we can get the condition number of the matrix sequence T J J * .The results are as follows in Table 2 (Rcond means the Reciprocal of the condition number).
Comparing the two results, we can see that before oscillation occurs, except 6 th -iteration has tiny difference, both results are exactly same.It shows, even if the condition number of the T J J * increases to the J's square, we can still get the good results.This fact proves again that bus-11 is not an ill-conditioned system.
What about the Point Diagonal Preconditioning E-Ψct for the gradient equation?
The differential equation should be: Here DD is the diagonal matrix of T J J * .We have the following results in Table 3.One may think 8858 gradient function evaluations are too expensive, but please note the step size is 20.0 -600.0.If we use fixed iteration (Euler method), the step size is not greater than 0.76.We try this method, the

Iteration
Step size Norm of gradient Norm of ( ) step size is from 0.73 − 0.76, it needs 12,902,751 evaluations.The proportion of the number of the evaluating gradient functions is 12902751/8858 ≈ 1456.In the later period of iteration, these two iterations keep monotone descending.When F reaches 0.0002089, LAVM reaches 0.0001027.They keep no variation.So we can de- termine the minimum of F is 0.0002089.At this time, can we say an efficient seeking minimum method, E-Ψtc, has been achieved?Some people can reject this statement, the reason is conventional Gauss-Newton method only needs 7 iterations, which can make LAVM reach 0.0002557.Ours answer is that, after 7 iterations, the LAVM continues to oscillate and if the system is large, say the dimension is 10,000, then to solve 7 linear equations with dimension of 10,000, the workload is surely much greater than 8858 function evaluations.

Discussion
Usually we evaluate the accuracy of a solution by the mismatch of the equation, but for ill-conditioned problem, it might be misleading.Here we give an extreme example.Consider a 2-dimension linear algebraic equations: If substitute x 1 by 0.9911, and x 2 by −0.4870, the mismatches will be −10 −8 and 10 −8 .People might think the most accurate solution has been obtained.Unfortunately, s/he is wrong.The exact solution of these equations is 2 and −2.From a mathematical theory, discriminating the accuracy of a linear equation, we need to resort the knowledge of the condition number [16], but this is a very difficult task to implement, especially for a real world problem.
We suggest taking the following easy ways to identify if the system is ill-conditioned or not: Using both single and double precision to do the calculation, if the results basically keep coincidence, we can consider the system is well-conditioned.
If using a Normal Newton method, for 43-bus system, single precision cannot make the LAVM less than 1E−5, but for double precision, the LAVM can go down less than 1D−12 easily, so 43-bus system is surely an ill-conditioned system.
For 11-bus system, we do the same process above, the results of single and double precision are almost the same, the LAVM cannot go down to 1E−4 and 1D−4.After 7 Newton iterations, the LAVM are 0.255709E−03 and 0.255653D−03, then oscillation occurs for both iterations.So we confirmed that the 11-bus system is not an ill-conditioned system, but the LAVM cannot go down to zero.
In order to comprehend this phenomenon, consider the equation: ( ) 2 0.0001 0 The solution does not exist, but at x = 0.0, ( ) In engineering design, making 0 F = or very small, e.g.BUT, finding the minimum of ( ) F x is an even harder task when comparing with solving ( ) 0 F x = .In the process of solving F x around some small value is oscillating, we can recognize that the minimum has been obtained.For 43-bus system, the minimum is about 1D−12, so this system is well designed.As for the 11-bus system, the E-Ψtc method gives the minimum 1D−04.We consider it is determined by the parameters of design, but it was said to be "operational" [1] pg. 1738.
To summarize our opinion, the 11-bus system is not ill-conditioned, the norm of ( ) F x does not go down to zero, the LAVM can reach 10 −4 , but in engineering it is "operational".For all the methods tested by us, the only method which can make LAVM to reach to 1D−04 is the E-Ψtc method.Now let us discuss the reliability of the algorithm.In many literature, there has been too much emphasis on measuring the efficiency of the algorithm proposed by the author(s) but not enough on reliability.For the nonlinear equations, the non-uniqueness of the solution makes measuring the reliability being an important subject.Reference [14] proposed a strategy to test the reliability: Giving a standard starting point x s which is regarded as very close to the solution, then choosing another two starting points 10x s and 100x s , which are regarded as being far away from the solution.This strategy is totally not working for the power system.Almost all people take flat start (1.0 + j*0) as the standard starting point x s , choosing (10.0 + j*0) and (100.0 + j*0) as the starting point is unimaginable.
For the power system, the authors of this paper propose the following strategy to test the reliability: 1) Using single and double precision to solve the equation ( ) 0 F x = .If the LAVM of both goes down to less than 10 −5 or 10 −6 and the results basically keep coincidence, we can believe the results are reliable.If the results do not coincident, like the 43-bus system, which reminds us that this system could be an ill-conditioned system.In order to get reliable results, we can use double precision to make even smaller LAVM.
2) If single and double precision are used, the iterations do not converge, the doubts about the algorithm and the existence of the solution might be raised.We suggest that people change the algorithm first to see what happens.If the iteration still does not converge, the existence problem should be considered.We observe all the iteration processes, if all of the iterations diverges, we can affirm that the solution does not exist.If the solution around some value oscillates and this value is relatively small, like the 11-bus system, we can recognize the solution exists, in engineering it is called "operational".We do not suggest readers to change the initial values to improve the convergence, because for the power system, flat start is good enough, these initial values are very close to the solution.
Even though oscillation can also be regarded as divergence, in this paper, we discriminate those two different phenomena.Here, so-called divergence means drastic oscillation happens in the iteration process.
3) As we have mentioned before, the solution of a nonlinear equation is not unique.The small LAVM does not guarantee the solution we searched is the solution what we want to get.The reliability asks us to list the results produced by different methods and to compare them, but none of [1] or [2] gives the final results.Especially [2] should more to this demand, because to test a new method, the reliability is more important than anything else.On the last page (pg.3657) of [2], there was a result of the 11-bus system, cited other literature.In order to compare, we list three results (voltage magnitude and angle): 1) page 3657, 2) optimal multiplier (6 iterations) and 3) point diagonal preconditioning E-Ψtc.
The value of the angle of the bus-2 is probably wrong.
From the three results above in Table 4, we can see they are basically coincident.The LAVM = 0.995493D−04 for E-Ψtc method, and the LAVM = 0.447367D−03 for 6 Optimal Multiplier iterations.As we pointed out before, after 9 optimal multiplier iterations, the LAVM keeps 0.2155D−03.In order to compare with the data on page 3657, we list the 6 iterations results instead of 9 iterations results.

Conclusions
In this paper, we give a review about references [1] and [2], and introduce E-Ψtc method as well as its application for Load Flow Problems.
We find out that the benefits of the optimal multiplier method can only occur in the initial period/stage of Newton iteration.For the 43-bus system, if the stopping condition is LAVM < 10 −8 , it can save one Newton iteration.For the 11-bus system, it can eliminate the oscillation which occurs in the normal Newton method.During the later period, we should take care of the loss of quadratic convergence.The applicability of the optimal multiplier method is only for the Newton iterative method.Other methods, such as steepest descent method, are not available.43-bus system is an ill-conditioned, using double precision can overcome "convergence trouble".11-bus system is not an ill-conditioned.It does not have convergence trouble.The problem is that the LAVM of ( ) F x is not as small as we usually expect, say 10 −5 or 10 −6 .It is likely 10 −4 .Here we use the word "likely" to express the idea we are not sure.This is because that we resort solving ( ) 0 F x = to get the LAVM.From Note: For the data on page 3657, in the front of angle should have a negative sign.
we can determine the right side vector b.The iterative initial make the Euclid norm less than 1D−08, it seems impossible.How about the Explicit Pseudo-Transient Continue" method.(E-Ψtc)?The following data can give the answer: (Variable step size)
x reaches the minimum 0.0001.Let us see what happens if the normal Newton method is used.The initial value of x 0 = 0.5, , after a sequence of monotone descend, at 6 th iteration, x = 0.2995D−02, ( ) F x = 0.1090D−03, then oscillation occurs, at 10 th iteration, x = −0.5698d−03,( ) F x = 0.1003D−03, then oscillation occurs again.This oscillation continues until ( ) F x reaches its minimum 0.0001, at this time the value of x is 0.8337D-04.If ( ) F x = x 2 = 0, the oscillation does not occur for Newton iteration.After 536 iterations, x = 0.2223D−161, ( ) F x = 0.4941D−323.After 537 iterations, x and ( ) F x keep 0.1111D−161, and 0.0000, have no any varia- tion.That is to say, we have an approximate solution, and the accuracy reaches 0.1111 × 10 −161 .

6 ,
usually the design is considered acceptable.Strictly speaking, for the engineering design problem, the solution of ( ) 0 F x = , in the sense of mathematics, does not exist.What we can do is to find the minimum.As for the size of the minimum, which is determined by the requirements of the engineering projects.

Table 1 .
In the table, NFE denotes Number of Function Evaluations.Double precision was used.
we use a fixed point iteration method to solve this linear equations The matrix B is diagonal dominant, but matrix A is not.This is why for every bus, we arrange the sequence according to ( ) , f e .Only by this way, we can guarantee every sub-block of the matrix is diagonal dominant.Here we would like to emphasize this matter again.The algebraic equations and differential equations were all the same, just the unknown variables have different order-one is ( ) located on the left half-plane.Choosing proper parameter, the fixed point iteration will converge.From the differential equations' point of view: For differential equations:

Table 2 .
Results of Newton Method for gradient equation by LINPACK DGECO.

Table 3 .
Results of Point Diagonal Preconditioning E-Ψtc for gradient equation.

Table 4 .
Comparison of three results.