Numerical Solution of Differential Equations | Power Systems - Electrical Engineering (EE) PDF Download

Numerical Solution of Differential Equations

Why numerical methods ? 
The reason we have taken a minor diversion here is that we wish to know the techniques to "understand" how systems described by differential equations behave. For very simple systems like the one below:

It is clear that x=0 is an "equilibrium" solution of the system, since, the LHS of the above equation equals zero at this value of x. In general, the behaviour of x, when its value at time t = 0 is x(0) is given by :

Verify that the above "solution" satisfies the differential equation. Note that the solution is in terms of a well known exponential function. If a > 0 , then the magnitude of x(t) keeps increasing with time if x(0) is not zero. On the other hand, if a<0, then x(t) tends to go to zero as time progresses. Thus, we are able to gain an insight into the behaviour from the solution given above.

However, it turns out that for many systems, it is not possible to write down the solution in terms of well understood simple functions. This occurs when the RHS of the differential equations have terms which are nonlinear or time variant functions of the variables. For example, the behaviour of rotor angle and speed deviation for a synchronous machine is described by the non-linear differential equations:

To understand the behaviour of such a system one has to turn to the numerical solution of these equations.

Numerical Solution of Differential Equations 
Since it is difficult to obtain the solution for a nonlinear set of differential equations, we try to utilize a computer program which numerically computes the solution at discrete points in time. An alternative would have been to implement a setup using scaled physical elements which mimic the differential equations, i.e., an analog computer. However, given their flexibility, numerical evaluation using computers is convenient and economical.

So how do we solve the differential equations numerically ? Let us consider a simple example:

If we wish to obtain the solution denoted by x0, x1, x2, ....xk at the closely spaced time instants t = 0, h, 2h, .....,kh, then we can approximate the above equation as:

Therefore if xk is known, xk+1 can be obtained by

Therefore, if the initial condition x0 is known, we can recursively find the solution of x at discrete instants of time. The approximation is likely to work well only if h is "small enough".

To illustrate this, consider the numerical solution at discrete points and compare it with the correct solution:

Let us take the constant a = -4.
We see from the figure above that the approximation works better for smaller h.
The method that we have presentated is in fact the simplest possible one and is called "Euler's Method" of numerical integration.

It is not difficult to see that another possible approximation of the differential equation could be:

This is called the "Backward Euler" method. The numerical solution obtained using different approximations and different values of h could vary. Inappropriate selection of the method and the time step h can lead to errors which can cause erroneous conclusions.

Comparison of the two methods s 
We have considered two methods: Euler Method and Backward Euler method. Euler's method is simple to implement on a computer program, but is known to have poorer numerical properties. In particular, it may wrongly show a system which is actually stable to be unstable. As an example consider the following equations:


If we use Euler Method, we descretise the equations as follows.

For a time step, h=0.1, with initial conditions x1(0)=1 and x2(0)=0.01, the numerically evalauated response is unstable, although this system is actually stable and has a time response given by:

One has to use a very small time step (say, h = 0.005s) if one has to avoid getting a grossly incorrect numerical solution. However, this would mean more computation time, even one may be interested only in capturing the slow response corresponding to the exp(-5*t) term in the response.



Let us repeat this exercise with Backward Euler with the same h=0.1. The discretised equations are:

Clearly, the solution obtained from Backward Euler Method does not "blow up". It also captures the slow transient corresponding to exp(-5*t) quite accurately with this value of time step.


The problem with Backward Euler Method is evident when we apply it to the solution of swing equations:

Clearly, getting the values of d and w at the k+1 th instant from the values at the k th instant is tough as we have to solve nonlinear algebraic equations (how does one solve them?).

Clearly there is a trade-off between numerical accuracy and complexity of implementation! We now consider a method which is not too difficult to implement and also gives reasonably good accuracy.

Runge-Kutta Fourth Order Method
Consider the differential equation:

The Runge-Kutta 4th order algorithm uses intermediate points in an interval of duration h to calculate the state xk+1 from xk as follows:

Note that the number of calculations per time step is larger than that for Euler Method. However, this method is more accurate than Euler method (for the same time step) for most systems.

Can you program an algorithm which implements this method for the swing equations ?