Numerical integration of nonlinear differential equations by use of rational approximatio

Fair, W.

Luke, Y. L.

Wimp, J.

English

NASA Technical Reports Server

I 
I 
i 
I 
' 1  
f 
Q 
4 
. ,  
NUMERICAL INTEGRATION OF NONLINEAR 
DIFFERENTIAL EQUATIONS BY USE OF 
RATIONAL APPROXIMATION 
FINAL REPORT 
15 January 1965 - 14 July 1967 
Task Order Contract No. NASr-63(07) 
NASA Hq. R&D 80X0108(64) 
MRI Project No. 2760-P 
GPO PRICE $ 
CFSTI PRICEtS) $ 
For 
Hard copy (HC) .7& 
Microfiche (MF) 
ff 663 July 66 
Office of Grants and Research Contracts 
Code SC 
National Aeronautics and Space Administration 
Washington, D. C. 20546 
11, -r f i  
425 VOLKER BOULEVARD/KANSAS CITY, MISSOURI 64110/AC 816 LO 1 -0202  
https://ntrs.nasa.gov/search.jsp?R=19670021846 2020-03-12T11:01:33+00:00Z
WRICAL INTEGRATION OF NONLINEAR 
DIFFERETVTUL EQUATIONS BY USE OF 
RATIONAL APPROXIMATION 
Yudell L. Luke 
Wpan Fai r  
Jet Wimp 
FINAL REPORT 
15 January 1965 - 14 July 1967 
Task Order Contract No.  NASr-63(07) 
NASA Hq. R&D 80X0108(64) 
MEU Project No. 2760-P 
For 
Office of Grants and Research Contracts 
Code SC 
National Aeronautics and Space Administration 
Washington, D. C.  20546 
(? jL-3 ” - -  c 1% c- a : .? 2s L- -3 I- 3 - 2  L3 
425 VOLKER BOULEVARD/KANSAS CITY, MISSOURI 64110/AC 816 L O  1-0202 
Aerona 
This report  covers research i n i t i a t e d  by Headquarters, National 
t t i cs  and Space Administration on Contract NASA Hq. R&D 80X0108(64) 
10-74-740-124-08-06-11, 3310-2487. 
Structures  .If 
complished a t  Midwest Research I n s t i t u t e  with m. Howard Wolko as pro jec t  
monitor. 
"Nonlinear EynCunics O f  Th,in She l l  
The research work upon which t h i s  report i s  based was ac- 
The authors take t h i s  opportunity t o  thank Miss RosexWry Moran 
for her ass is tance.  
Approved for: 
MIDWEST RESEARCH INSTITUTE 
Sheldon L. Levy, D i r  
Mathematics and %ys 
14 July 1967 
. 
l i d  
TABLE OF CONTENTS 
Page No. 
. . . . . . . . . . . . . . . . . .  I. Summary and Introduction 1 
11. Development of Integration Formulas . . . . . . . . . . . .  2 
111. Truncation Error . . . . . . . . . . . . . . . . . . . . . .  5 
IV. Examples and Applications . . . . . . . . . . . . . . . . .  8 
v. Conclusions . . . . . . . . . . . . . . . . . . . . . . . .  10 
References . . . . . . . . . . . . . . . . . . . . . . . . . . . .  12 
iii 
S. SUMMARY AND INTRODUCTION 
equati 
In the 
ns there 
numerical treatment of ordinary linear differential 
re many effective integration formulas which yield quick, 
reliable results. These techniques depend heavily on the fact that solu- 
tions to linear differential equations are very well behaved in that all the 
needed higher order derivatives exist in the range of interest. 
In general, the theory of ordinary nonlinear differential equa- 
tions offers no clue as to the singularities of the solutions of such 
equations. 
heuristically. Obviously the usual numerical integration techniques fail 
in the region of such singularities. Not only are the results inaccurate 
near a singularity, but also the location of such a point evades detection. 
Hence, new techniques must be developed which will deal effectively with 
the problem of singularities of solutions to nonlinear differential 
equations. 
Thus the detection of singularities must be accomplished 
Recently we obtained an algorithm for computing rational ap- 
proximations to the solution of a wide class of nonlinear differential 
equations, see [l] . Although this technique cannot be extended to arbi- 
trary nonlinear equations, it clearly demonstrates the power of rational 
approximations in dealing with a function which has singular points in the 
range of interest. These approximations ascertain the existence of zeros 
and poles of a function and locate these critical points with great ac- 
curacy. Rational approximations also allow accurate computation of the 
function near a singular point -- a decided advantage over the usual ap- 
proximations. 
In light of the advantages cited above, it is desirable to obtain 
formulas for numerical integration of ordinary nonlinear differential 
equations which are based on rational approximations. 
In Section I1 we develop multi-step predictor and corrector 
formulas for numerical integration. In Section I11 we obtain the truncation 
errors associated with such formulas. The results are applied to some ex- 
amples in Section IV. A recommended procedure for the use of these rational 
approximations in the numerical treatment of nonlinear differential equations 
is outlined in Section V. 
11. DEVELOPMENT OF INTEGRATION FORMlTLAS 
Here we develop multi-step predictor and corrector formulas based 
on rational functions for the numerical integration of the differential 
e quat ion 
The extension of the formulas to higher order equations is immediate. 
It is assumed that all the needed starting values are available. 
We wish to approximate the solution to (2.1) by a rational function. 
write 
We 
where the coefficients aj and bj are unknown. 
In order to obtain a k-step method, it is sufficient to require 
that either n+m = 2k , or n+m = 2k-1 . For convenience, we set 
and y(xj) = yj 
We develop the predictor and corrector formulas for the case m+n = 2k and 
then indicate the necessary changes for the case 
for a set of equally spaced interpolation points 
xo<xl<...<xk+l . Without l o s s  of generality, xo = 0 and xj+l = Xj+l . 
m+n = 2k-1 . 
Let m+n = 2k and require that in (2.3) - (2.4), 
2 
h,n(xj) = 0 j = 0,1, ..., k+l , 
The resulting system of equations in the unknowns aj and bj has a 
solution, if and only if, the following determinant vanishes; 
Ap = 
Rk+l,m+l 'k+l, n+l 
Tk,m+l 'k , n+l 
= o ,  
where the entries in the determinant are rectangular arrays of the size 
rectangular array 
indicated by the subscripts, and if hij denotes the i,j- th entry in a 
Kp,j , the arrays in the determinant (2.6) are defined 
by 
(i-l)j-', i = 1,2, ..., k+2 ; j = 1,2, ..., m+l and rll=l , rij = 
sij = (i-l)j-%i-l, i = 1,2,.. .,k+2 ; j = 1,2,. ..,n+l and s l l ~ o  , 
(2.7) 
tij = (j-l)(i)j-*, i = 1,2 ,..., k ; j = 1,2 ,..., m+l , 
u )= (j-l)(i)3-2yi+(i) j-1 yi, i = 1,2 ,..., k ; j = 1,2 ,... ,n+l . 
ij 
We call Ap the predictor determinant. To obtain the corrector determinant, 
A c  , corresponding to the conditions (2.5), replace the first row of (2.6) 
by the last row of (2.6) with k replaced by k+l  . After setting 
ALP = Ac = 0 and expanding the determinants, we obtain 
3 
where the  coe f f i c i en t s  i n  (2.7) depend on y j  
and the  coe f f i c i en t s  i n  (2 .8)  depend not only on these  values,  bu t  a l so  on 
y'k+l . 
corrected value. This i s  discussed i n  more d e t a i l  l a t e r .  Thus a f t e r  s e t -  
t i n g  
rec ted  values a t  yk+l by solving (2.7)  and (2 .8) .  The predicted value of 
yk+l  
a predicted 
s t a b i l i z e .  However, f o r  small h , t h e  pred ic tor  can be used alone. 
and y l j  , j = O,l, ..., k , 
Note t h a t  (2 .8)  requires  solut ions of a quadrat ic  equation f o r  t h e  
y ( x )  = y(xo+hxj we obtain the f i n a l  Torm f o r  t he  predicted and cor- 
i s  obtained from (2 .7) ,  t h e  d i f f e r e n t i a l  equation i s  employed t o  ge t  
ylk+l and the  corrector  i s  used repeatedly u n t i l  t h e  r e s u l t s  
For t h e  case m+n = 2k-1 , we ge t  t he  pred ic tor  from ( 2 . 6 )  by 
de le t ing  the  f i r s t  row and we g e t  t h e  corrector  from t h i s  r e s u l t  by in- 
creasing t h e  row indices  of  the  l a s t  k rows by uni ty .  In  i l l u s t r a t i o n ,  
t h ree  cases a r e  l i s t e d  below. 
If m=n=l : 
2 
2yoYl -5 l+~YoY~ '  
S 0 - 5 1 + h Y 1  1 '  
y2(predic ted)  = 
and t h e  corresponding cor rec tor  formula i s  
2 2 
Y2 -2y1Y2 + ( Y 1  4Y11Y2') = 0 
For m = l  , n=2 : 
and f o r  t he  cor rec tor ,  t he  values i n  (2 .8)  a r e  given by 
4 
For m=n=2 : 
where 
d3 = hY ' (4Y -3hY '-4Y2) - (Y1-Y2)(5Y1-5Yo+6hYo' ) Y 1 0  0 
and for the corrector, the values in (2.8) are 
1 c = Y1-Yo-hYo'- 4 h Y l '  , 
D = 2hy0'y1+yo2-yl2 +Yl-YO+ 5 hYOY1' 1 9 
2 1 
E = YoY12-Yo Y1-hYO'yl2-hy2' (Y1;Yo)2 - - 4 hyo?yl'+h3y0'y1'yp' . 
111. TRUNCATION ERROR 
In this section we develop an expression for the truncation error 
associated with the integration formulas of Section 11. We restrict our- 
selves to the predictor formula for the case m+n = 2k (see (2.5) - (2.7)). 
The other cases are handled similarly and, since we are interested only in 
the order of the truncation error, the treatment of one case is sufficient. 
As in (2.2), we set 
Let x,<x <...<x~+~ be equally spaced interpolation points and assume 
(2k+l) times differentiable in the interval [ ~ ~ ~ x ~ + ~  ] except at the poles 
C1,C 27...,Cv of orders s1,s2, ..., sv respectively. Thus, if 
y is 
5 
J 
j=1 
Then r ( x ) y ( x )  i s  (2k+l) times d i f f e ren t i ab le  i n  [ x ~ , x ~ + ~ ] .  Now we s e t  
where 
k 
j =1 
n ( x )  = TT (x-xj)  2 , 
m, n 
and F(x)  i s  a f'unction t o  be determined. Now i f  
h , n ( X )  = Y(x)$(x) - Pm(x) 9 
then f o r  t he  case m+n = 2k , t h e  predictor  formula requires  
&,n(xi)  = 0 , i = 0,1, ..., k+l , 
& , n ' ( ~ i )  = 0 , i = 1,2 ,..., k . 
Consider t he  funct ion 
We have 
g(xi) = 0 = g ' (x i )  , i = 1 ,2  ,..., k . 
6 
(3.4) 
(3.5) 
(3.7) 
Now c and d i n  (3 .6)  can be determined so  t h a t  g ' (xo )  = 0 . Thus, 
g " ( t )  has 2k zeros i n  the  in t e rva l  [Xo,xk+l] . Repeated use of Rol le 's  
theorem guarantees t h e  exis tence o f  a number 
x ~ < U < X ~ + ~  such t h a t  g(2k+1)(u) = 0 . Differen t ia t ion  of (3 .6)  (2k+l)  
times y i e lds  
u i n  the  i n t e r v a l  
Solving f o r  F (x )  , we ge t  
so t h a t  
(3.10) 
Using the  f a c t  t h a t  x = Xo + j h  , s e t t i n g  x = x k+ l  and simplifying, we 
have j 
(3.11) 
Thus, t he  order of t he  t runcat ion e r r o r  i s  
r ec to r  i s  t h e  same. It i s  very important t o  point  out t h a t  t h e  usual  
l i n e a r  k-step in tegra t ion  methods have t runca t ion  e r r o r  of a t  most order 
of  Thus, the  predictor  (2 .7 )  i s  more accurate than t h e  
usual  k-step method. 
h2k . The order f o r  the  cor- 
hk+* , see [2] . 
7 
For the  case n+m = 2k-1 , it can be shown as above t h a t  t he  
predictor  and corrector  both have truncation e r rors  of order 
f o r  computational purposes, it seems desirable  t o  u t i l i z e  the approximations 
(2.6) i n  which 
h2k . !Thus 
n+m = 2k-1 . 
Iv. EXAMPLES AND APPLICATIONS 
We apply the  r e s u l t s  of Section I1 t o  three  examples. In the  
f i r s t  two examples, f i r s t  order d i f f e ren t i a l  equations a re  t r ea t ed .  Here 
t h e  approximations (2.7) - (2.8) w i t h  m = l  and n=2 are  employed. I n  the 
t h i r d  example, a second order d i f f e r e n t i a l  equation, we use m=n=2 . In  
both cases, the corrector  i s  used repeatedly u n t i l  successive i t e r a t e s  
agree t o  e ight  decimal places.  
bes t  w i t h  the  previously computed value of 
The root  of (2 .8)  was chosen which agreed 
yk+l . 
Let u = tan(x+n/4) , v = Jl(x)/J0(x) and l e t  Z be Painlev6's 
second transcendent. These functions s a t i s f y  
1 v '  = l+v2 - - v , v(0) = 0 , 
X 
and 
Now u , v and Z have simple poles a t  n/4 , 2.40482 and 1.1577, 
respect ively.  See Tables I, I1 and I11 fo r  the r e s u l t s  of numerical in te -  
grat ion.  
t h a t  these approximations do indeed detect  the presence of a pole i n  each 
case. To i l l u s t r a t e  t h i s  l a s t  remark, we have l i s t e d  i n  Table I1 t he  ap- 
proximations t o  the smallest pole of tan(x+rr/4) based on use of t h e  
predictor  only. 
quadratic,  and the  zeros of the  denominator were computed. We point out 
t h a t  t h i s  t ab le  a l so  indicates  t ha t  t he  pole can be computed very accurately 
i f  the  funct ional  values a re  known r e l a t i v e l y  close t o  t he  pole. 
Notice, t h a t  the r e l a t ive  e r ro r  has a slow r a t e  of growth and 
Here the  function w a s  approximated by a l i nea r  over a 
I n  Tables I, I11 and  I V ,  t he  corrected value t h a t  appears i s  t h e  
f i n a l  r e s u l t  of correcting repeatedly u n t i l  successive corrected values 
agreed t o  e ight  decimals. For the s tep  s ize  h = 0.01 , a l l  t he  interme- 
diate values were computed but these a re  omitted f o r  t he  sake of brevi ty .  
Note t h a t  the  computed values beyond the  pole a re  very accurate. 
though t h i s  pract ice  of integrating over a pole i s  not recommended, it does 
give valuable information about the behavior of the  function: 
Even 
TABLE I 
NUMERICAL SOLUTION O F  U' = l+u2, u(0) = 1 
u = tan(x+rr/4) 
h = 0.05 h = 0.01 
- X True. (u) Predicted Corrected Predicted Corrected 
0 .1  
0.2 
0.3 
0.4 
0.5 
0.6 
0.7 
0.8 
0.9 
1 .0  
X - 
1.22305 
1.50850 
1.89577 
2.46496 
3.40822 
5.33186 
11.68137 
-68.47967 
- 8.68763 
- 4.58804 
1.22304 
1.50848 
1.89574 
2.46493 
3.40815 
5.33165 
11.67998 
-68.5 966 7 
- 8.73393 
- 4.62137 
1.22305 
1.50850 
1.89577 
2.46498 
3.40826 
5.33195 
11.68153 
-68.66273 
- 8.68629 
- 4.64804 
1.22305 
1.50850 
1.89577 
2.46496 
3.40822 
5.33186 
11.68138 
-68.48685 
- 8.69860 
- 4.56120 
1.22305 
1.50850 
1.89577 
2.46496 
3.40822 
5.33186 
11.68139 
-68.49443 
- '  8.69493 
- 4.56121 
TABLE I1 
LOCATION OF POLE OF tan(xcrr/4); TRUE POLE; rr/4 = 0.785398 
h = 0.05 
Root of  Q ~ ( x )  
0.60 
0.65 
0.70 
0.75 
0.80 
0.7868 
0.7855 
0.7869 
0.7851 
0.7851 
X - 
0.76 
0.77 
0.78 
0.79 
0.80 
h = 0.01 
Root of Q ~ ( x )  
0.78540 35 
0.78539 82 
0.78540 08 
0.78539 89 
0.78539 65 
TABLE I11 
X - 
0.2 
0.4 
0.6 
0.8 
1 .0  
1.2 
1 .4  
1 .6  
1.8 
2.0 
2.2 
2.4 
2.5 
1 
ST v, v(0) = 0 NUMERICAL SOLUTION OF v' = l+v2 - 
v = JI(x)/J0(x) 
h = 0.05 h = 0.01 
Corrected True ( v )  Predicted Corrected Predicted 
0.10050 
0.20411 
0.31436 
0.43584 
0.5 7508 
0.74246 
0.95606 
1.25141 
1.71041 
2.57592 
5.03762 
207.43659 
-10.27398 
0.10050 
0.20411 
0.31436 
0.43584 
0.57508 
0.74246 
0.95606 
1.25141 
1.71041 
2.57589 
5.03636 
126.9513 8 
-11.13360 
0.10050 
0.20411 
0.31436 
0.43584 
0.57508 
0.74246 
0.95606 
1.25142 
1.71043 
2.57596 
5.03698 
128.87352 
-11.64502 
0.10050 
0.20411 
0.31436 
0.43584 
0.57508 
0.74246 
0.95606 
1.25141 
1.71041 
2.57592 
5.03754 
207.31622 
-10.28250 
0.10050 
0.20411 
0.31436 
0.43584 
0.5 7508 
0.74246 
0.95606 
1.25 141 
1.71041 
2.57592 
5.03754 
207.31021 
-10.28618 
TABLE IV 
NUMFRICAL SOLUTION Z" = 2Z3+xZ+1, Z(0) = 1, 2' ( 0 )  = 0 
Z IS PAINLEVE'S SECOND TRANSCENDENT 
h - 0.05 h = 0.01 
- X !rrue (z) Predicted Corrected Predicted Corrected 
0.2 1.06261 
0.3 1.14638 
0.4 1.27415 
0.5 1.45921 
0.6 1.72538 
0.7 2.11844 
0.8 2.73694 
0.9 3.83440 
1 .0  6.31100 
1.1 17.31546 
1 . 2  -23.64085 
1.06271 
1.14634 
1.27377 
1.45 751 
1.72170 
2.11211 
2.72608 
3.81512 
6.24525 
19.30070 
1.62902 
1.06267 
1.14640 
1.27379 
1.45730 
1.72159 
2.11185 
2.72581 
3.81417 
6.25787 
21.69210 
-3.06227 
1.06261 
1.14637 
1.27415 
1.45921 
1.72537 
2.11840 
2.73710 
3.83522 
6.31763 
1 7 . 3  7845 
* 
1.06261 
1.14638 
1.27415 
1.45921 
1.72538 
2.11845 
2.73710 
3.83520 
6.31758 
17.38471 
* This value could not be computed due t o  overflow on t h e  computer. The 
las t  value computed w a s  f o r  x = 1.14. 
10 
V. CONCLUSIONS 
It is evident that the integration techniques developed here can 
he extremely useful, particularly in locating singularities of solutions to 
nonlinear differential equations and in computing functional values near 
such points. Before these approximations can be used on a wide scale, 
several problems must be analyzed. First there is the problem of selecting 
the correct root of the quadratic equation (2.8) associated with the cor- 
rector. Another pertinent question is the stability of the process. This 
is particularly important when the predictor is used without the corrector. 
Finally, there is the problem of obtaining stepvise a priori estimates of 
the error involved in the integration. In connection with the latter, it 
appears that the predicted value and the first corrected value can be used 
to estimate the error in much the same way as is done f o r  :he usual 
numerical integration schemes. In the present instance the process is 
much more complicated for at each stage one must compute the coefficients 
in the denominator polynomial of the rational approximation and other 
quantities as evidenced by (3.11). 
and we recommend further research on these topics. 
A l l  these questions are interrelated 
Our experience to date with the developments in this report sug- 
gests that the approximations be employed in the following fashion. 
usual integration formulas should be used as long as the solution to a 
given differential equation is well behaved. If large differences in suc- 
cessive computed values occur (indicative of a pole), employ the rational 
approximations to locate the pole p . Then introduce a reciprocaltrans- 
formation in the differential equation to obtain the equation whose solution 
has a zero at p . In the neighborhood of this zero, the usual integration 
techniques can be utilized. After having gone a sufficient distance beyond 
this zero, use the reciprocal transformation to recover the original 
equation and proceed with the integration. 
The 
11 
REFERENCES 
1. Fair, W., and Y. L. Luke, "Rational Approximations to the Solution of 
the Second Order Riccati Equation," Math. Comp., Vol. 20, No. 96 
(1966). 
2. Henrici, P., Discrete Variable Methiods in Ordinary Differential 
Equations, Wiley and Sons (1962). 
. 
. 
12 


Numerical integration of nonlinear differential equations by use of rational approximation Final report, 15 Jan. 1965 - 14 Jul. 1967

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server