Thayer Watkins Silicon Valley & Tornado Alley USA |
---|
Centered Difference Approximations of Derivatives in the Numerical Solution of Differential Equations |
The centered difference approximation to the first derivative has the advantage of greater accuracy over the forward and the backward difference approximations but it necessarily introduces an additional solution, called the computational mode, which does not correspond to a solution of the differential equation. If this computational mode grows in magnitude as the numerical solution is computed it may make the numerical solution worthless as an approximation to the solution of the differential equation. The analysis below points out that the growth of the computational mode requires an initial deviation from the true solution and this makes the computation at the intial steps of the solution crucial. Since the centered difference introduces a higher order for the difference equation than that of the differential equation its solution requires more initial conditions and those additional initial conditions have to be computed using some approximation other than the centered difference. The errors introduced by those initial condition approximation may grow and propagate to the detriment of the accuracy of the later computations.
Consider the simple differential equation
This is a linear first order differential equation and has the solution
The centered difference approximation is
where h is the time step Δt and y_{i}=y(ih). This is a second order difference equation. Since it is linear and homogenous it has two solutions of the form y_{i} = λ^{i}. The parameter λ must satisfy the condition
The two solutions to this quadratic equation are
When h is sufficiently small such that (1+(hc)^{2})^{1/2}=1+(hc)^{2}/2 these solutions are closely approximated by
If
When z is sufficiently small ln(1+z) is closely approximated by z. Thus one solution corresponds to
In the limit as h goes to zero this goes to c and corresponds to the true solution to the differential equation. The other solution λ = -1 + hc(1-hc/2) is a cyclic solution of period 2. When c is positive the magnitude of this solution is less than unity and the difference scheme would be stable. If c is negative the magnitude is greater than unity and the solution would blow up; i.e., the deviation from the true solution would growl exponentially.
Now let us consider the process for getting the solution to the difference scheme. Although the numerical values are computed iteratively the process is equivalent to fitting the general solution of the difference equation to the initial conditions. The general solution is
where q_{1} and q_{2} are constants to be determined from the initial conditions; i.e.,
Let λ_{1} denote the physical mode and λ_{2} the computational mode. The computational mode is no problem if q_{2}=0.
While y_{0} is known and equal to the value of the true solution at time equal to zero, y(0), the value of y_{1} is not known because the true solution at t=h, y(h), is generally not known. Thus y_{1} usually estimated by some method outside of the iteration scheme.
The values of q_{1} and q_{2} are found as the solution to the equations:
The solutions are
From this we see that q_{2} will be zero if
Thus the ideal value for y_{1} is not the true solution value at t=h, y_{0}e^{hc}, but instead y_{0}λ_{1}. The deviation from the physical mode is directly proportional to the deviation of the value used for y_{1} from y_{0}λ_{1}. To the extent that λ_{1} is a good approximation of e^{hc} the true value of the solution at t=h is a good approximation of the ideal value of y_{0}λ_{1}.
The points above can be nicely illustrated. First consider the case of c=1. This is the stable case. The numeric solutions based upon centered difference but using two slightly different estimates of y_{1} for h = 0.1 and for 10 time steps are so close it is difficult to distinguish one from the other and both lie on the line of the true solution. This would also hold true if the number of time steps were increased to 100.
For the unstable case of c=-1 the match up of the centered difference solution for the two estimates of y_{1} (one computed using a forward difference and the other computed using a backward difference) is quite close to each other and the true solution for the first 10 time steps, as shown below.
But when the number of time steps is increased to 50 the picture is quite different. The numberic solutions begin to oscillate above and below the true solution to a drastic extent.
But the numeric solution based upon y_{1} being estimated by by y_{0}λ_{1} is right on the true solution throughout the 50 time steps.
When the number of time steps is increased to 100 the numeric solutions based upon centered differences and y_{1} from forward and backward differences blow up and run off the scale of the graph. But the numeric solution based upon y_{1}=y_{0}λ_{1} is right on target through the 100 steps, as shown below.
A note of caution: If the value of y_{1} is not exactly equal to the exact value of λ_{1} times the exact value of y_{0} the solution will eventually blow up no matter how impressively close it is during the initial time steps.
HOME PAGE OF Thayer Watkins |