CONTRIBUTION TO THE THEORY OF NUMERICAL
INTEGRATION OF NON−LINEAR
DIFFERENTIAL EQuATIoNs（III）
BYTATsuJIRo SHIMIZU
1 11．Error fbr pre面ctor雫oπ㏄tor method． When the local tnmCation error i皿一 creases hl the course of numerical integration of ordinary differential equations， we ◆ften continue to calculate the solution by decreas血g the step sセe． Henrici discusses m 1亘s book：（Discrete variable methods in ordinary differential euations pp．256−258）a mathematical ground to estimate the truncation error for some Unear multi−step methods（the predictor−corrector methods）． His book contains only two cases． The first case：The mlmerical valuesみcoi皿一 dde with the values of the true ’ 唐盾撃浮狽奄盾氏@7’（xn）（this assumption is not practical）at JC＝Xl， X2， ・・．， X”， °… The s㏄ond case： α為一1＊＝αカー1，α克一2＊＝α鳶．2，…，αo＊＝αo， whereα＊are the coef五一 くients ofγof the predictor， and αare the ’ coeMcients of y of the corrector． For the above two cases he shows that ア。＋ry。＋誤＝〃＋ユ（Cp．、＊−Cp．、）ア（坪）（x。）＋0（hP’2）． Adams−Bashfbrth’s method belongs to the second case， but the血mous Mine’s Emethod does not． We shaU show the above relation here fbr the general case． Let the di丘brential equation be（35） 農一∫（・，・）
where∫（x， y）satis丘es properly given assurhptions． Let the predictor be （36） γヵ＋為＊＋α鳶＿・＊γn＋ヵ＿・＋…＋α。字yヵ一乃｛β鳶r・＊γノ轟＋レ・＋…＋β・＊」〆ヵ｝ 寸he corrector be （37） yヵ＋k＋α鳶＿1アヵ＋ヵ＿、＋… ＋αo yn ’＝h｛β乃ア’ヵ＋k＋β為＿1y’抱＋h＿、＋… ＋βoyヵ｝ where both（36）and（37）have the local truncation e宜or of the order乃力＋1． Let the true solution of（35）be y（x）， then we have （38） 【39） y（x十kh）＋侮．、＊ア（x＋ん一1乃）＋…＋α。＊y（x） 一乃｛βh＿1＊γ’（x＋k−lh）＋… ＋β。＊γ’（x）｝＋Cp＋1＊ア（♪＋1｝（x）〃＋1＋0（乃ρ＋2） y（x十kh）＋α乃．、y（x＋丘一1乃）＋…＋α。γ（x） ＝乃｛β口’（x＋kh）＋β鳶＿、ア’（x＋丘一1乃）＋… ＋β。γ’（x｝＋Cp＋1夕（♪＋1｝（x）ht＋1 ＋0（hP＋2）． ［51］52
T．SHIMIZU
Now consider the fbllowing lhlear multi−step Inethod． （40） γお袖＊＊十ab＿1γヵ＋k＿1十・∴十αoγπ＝乃｛β彦＿1＊ア’n＋ヵ＿1十…十βo＊γヵ’｝ then ． （41） y（x十kh）十αk＿ly（x十瓦一lh）十…十αoγ（x） 一乃｛βfi＿、＊〆（x＋瓦一lh）＋… ＋β。＊γ’（x）｝＋Cク＋、＊＊ア（P＋1）（x）hP＋1＋0（乃♪＋2） We assume that en＝アーγ（xn）＝hPe（xn）十〇（hP＋1）， p≧1， after He面ci， where e（x）is acontinuously differentiable fUnction of x． Now en being the solution of a丘nite dfference equation， which is used to solve 加merically the differential equation fbrθα）we have θ。−hPe（Xn）＋h担q（x。）＋0（乃力＋2） ’ under the asspmption that∫（x，ア）and the h丘tial values are all analytic and the or− der of accuracy of the口Ut輌al value is larger than P十1． ・ Now we have consider the relations（42），（43），（44）which are obtained by（36）一 （38），（37）一（39），（40）一（41）respectively． γo，γ1，…，γぷ一1being all the same fbr（42），（43），（44）， onlyア勃＋k，γヵ＋為＊，γ勃＋ヵ＊＊ may be different each other． Now γ。．i’一〆（x。＋i）一∫（x。＋i，γ。＋i）一∫（x。＋i，ア（Xn＋i））−hP9（x糾i）θ。＋i，where
9（Xn＋i）＝fy（Xn＋i，γπ＋i） and eヵ＝hP（e（x．）十〇（乃））． In the fbllowing 3 relations（42），（43），（44）， em， e（xm）， gαが）are the same fbロη＝・ 0，1，2，…，n十ん一1． But en＋k， e（Xn＋k）， g（xn＋k）may be di丘’erent． Shlce these eヵ＋k， e（x潟＋k），9（xヵ＋ヵ）apPear only in the right hand side of the rela− tion （43）， the rdations （42），（43），（44） can be written as， （42） （43） （44） θπ＋k＊一一（αh＿、＊en＋ヵ＿、＋…＋α。＊εヵ）＋乃｛β，＿、＊9（Xn＋ヵ＿、）en＋為＿、 ＋…＋βo＊9（Xn）θヵ｝＋Cρ＋、＊ア（P＋1）（xπ）hカ＋1＋0（hP＋2） e，、＋鳶一一（αk＿1θち＋み＿1十…十aoeπ）＋h｛βk9（X針k）en＋k＋βゐ＿19（Xπ＋為，．、）θ糾為＿1 ＋… ＋β09（Xn）θπ｝＋Cρ＋、γ（カ＋1）（x鰺）hP＋1＋0（h？＋2） en＋カ＊＊一一（α為＿・e匁＋k＿・＋…＋α・en）＋乃｛β鳶＿、＊9（僻ヵ＿、）en＋為＿、＋…＋β。＊9（Xヵ）e＃｝ ＋Cρ．、＊＊γ（カ＋1）（Xn）hカ＋1＋0（hP＋2）． Now by 9（Xn＋i）en＋i−9（x。＋，）hP（e（x。＋．；）＋0（旬）−hP9（x。）e（x。）＋0（h餌） i＝1，2， … ， k−1，k． （42），（43），（44）gives respectively． （45） （46） （47） い Since these multi−step methods satisfy the convergent θ露＋為＊一一（偽＿1＊en＋カ＿1十…十α0＊en）＋乃力÷1｛（β必＿1＊＋…＋β。＊）9（Xヵ）θ（X勃）｝ ． ＋（：P．、＊y（P＋1）（x。）乃担＋0（乃力＋2） er．＋克一一（αk＿1θ％＋k＿1十…十αoen）＋乃力＋1｛（β左＋β丘＿、＋…＋β。）9（Xn）e（Xヵ）｝ ＋Cρ．、γ（ρ＋1）（XPt）hP＋1＋0（乃力＋2） en＋k＊＊一一（ak＿・＊en＋k＿、＋… ＋α0＊en）＋hP＋1｛（β，一、＊＋… ＋．β。＊）9（X”）θ（X．）｝ ＋Cp．、＊＊γ（餌）（Xn）hP＋1＋0（hク＋2）． condition of Dahlquist weCONTR田UTION TO THE THEORY OF NUMERICAL INTEGRATION、フ．．二53
have
k＋k−1ak−・＋…＋α・一β鳶＿・＊＋…＋β。＊一β鳶＋β，＿、＋…＋β。． Hence we have by（46）一（47） （48） e元＋鳶一θヵ＋鳶『c＊＝（Cp＋1−（］p＋1＊＊）γ（♪＋1｝（xヵ）乃ク＋1十〇（hρ＋2）． On the other hand， we have α‘θ頑一α‘（hPe（x。＋言）＋乃ρ＋1q（x。＋∂＋0（〃＋2）） ’ 一α∠（”e（x。）＋ihカ＋1e’（x。）＋0（〃＋2）＋hP＋’q（x。）＋0（hP＋2））． By the condition of convergence we have 1＋ak−・＋…＋α・−1＋αk．、＊＋…＋α。＊−0， え＋瓦一1α鳶．、＋…＋a、一ん＋ん一1αカ．、＊＋…＋α、＊． Hence we have by（47）一（45） （49） en＋牽＊＊−eヵ＋鳶＊＝（（7ρ＋1＊＊−Cp＋1＊）ア（ク＋1）（xヵ）ノip＋1十〇（乃ρ＋2）． Hence by．（48）一（49） （50） θ…一θ・＋・＊一γ。＋・一γ。＋・＊一（（7P．、＊一（］P．、）乃力＋ウ（カ・・｝（Xn）＋0（hP・・）． which is the relation required． “ This relation（50）is of theoretical one， since the step size is not so small in prac・ tice and O（hρ＋2）is not necessarily ・ small compared to the other terms． 12．D逝culty in numerical method． Runge−Kutta，s method or Taylor’s expan− sion method is used successfUlly to solve the血tial value problems of ordinary differential equations in some cases． But we meet o㏄asionally with some initial value problems， which can not be obtahled sat’isfactorily by any numerical methods． 丘1this s㏄tion we shall discuss the latter problems． There are two different stand points fbr numerical integration of ordinary difi7er− ential equations． The first case is the fbllowing：Assuming the eXistence of the solution on O≦x≦X， it is to integrate it numerically up to X from O． 、 シ ． This case has been studied by many writers． The truncation errors are estimated（though those are over−estimated）， and if we can take the step size of the numerichl mtegration so sma皿， then the numerical solu− tion y％at x＝xn，0≦x，、≦X can be obtained so that l yn一γ（xn）1＜ε fbr an arbitrariy smal1ε， whereア（xn）is the value of the true solution of the differential equation at x＝xn（without considering round off errors）． 、 If we can increase the number of d㏄imal digits as we please， then we can obtain γ外so that l yn一γ（xn）1＜ε， even under the consideration of round ofF errors． （C￡T．Shimizu：Co皿tribution to the theory of numerical integration of ord血ary differential equations， II， TRU Mathematics， Vol．4，1968）． B・f・・ew・di・cuss these p・・bl・m・，，we a・e tO・xp1・祖・・m・i・・tabmty・f n㎜・ricaI solution． Fox hl his book：（Numerical solutions of ordinary and partial differential equa一句 54
T．SHIMIZU
tions）discusses the equationめ7／dx＝アーx， and shows that the solution y＝x十1can not be obtailled．by numerical integration． He caHs such an instabMty‘‘㎞herent instabiHty・” 、 Further he discusses the e（luationばγ／dx＝−1qy十9−10x（the true solution bei皿g y・＝Ae−10x−x十1）． If the step size h＞0．2， the numerical solution fbr some Taylor’s expansion method increases without bound． He calls such an instability‘‘partial instabVUty．”39） It Inay happen that if the step size is not su伍ciently small a numerical method gives another solution which behaves quite differently from the required one． We shall call here such an instab伍ty‘‘partial instab伍ty”also． This case occurs when the required solution is inherently unstable or behaves sirnilarly． Some equation may have a solution which increases rapidly in the neighbourhood of the initial point x＝＝O． For example， fbrめ’／dx＝ア2 the solution y＝＝−100／（1十IOOx）withγ＝−100 atκ＝O can not be followed by numerical integration unless the step size may be taken so small（depending on the initial value）． If we take乃＝1／10 and apply Euler’s method， thenア（h）＞0， which is not the value of the above solution at x＝h． This corresponds to partial instability． The s㏄ond case： In general the solution of differential equations is not known to exist how侮r from x＝0． We are usuaUy to calculate the numerical solutionγカ by choos血g properly the step size hn． （This case occurs when we are study㎞g the non−linear osc皿ation theory and other important problems）． To solve numericaUy the differential equation， whose sohltion is not known， belongs to this case． ’For some dffereptial equation the solution may not or may have some movable singularity depending on the initial value． For example the equation tか／dx＝y2， whose祖itial value y−Yo at x＝0， has the solution y＝γo／（1−xyo）， which has a pole at x＝1∫γo． Even when the solution exists fbr O≦x＜oo， the truncation errors（also round off errors）increase in the course of numerical integrati皿fbr any numerical method （stable or not）． 39） Fox’s hlstability occurs while乃is not so large． But the equatiol1〆＝−1qγ十90−10x has a Iarge coe丘icient ofア． Letγ’＝−Ly withア（0）＝アo． If hL≧3， then Euler’s m・・h・d・・T・yl・d’・m・・h・d・・＋・一・醐’・芸・・〃・…・芸〆〃・・…th・・am・…ta− biHty wh且・h≦100 if L≧100・ Thus however small h may be and however high order of Taylor’s method is adopted， there may exist some differential equationア’＝−Ly， showing Fox’s partial instability． …T・y1…exp・n…nb・・・…＋・一・・（ （hLl−hL十 、 2）2＋…・・−1）・（響）一・…w・ can take hL so that IgKl＞1 fbr any given K， and L so large that h is smaller than any given Ilulliber． For such K and h Fox’s instability㏄curs．CONTRIBUTION TO THE THEORY OF NUMERICAL INTEGRATION
55 ノ Inherent instabi且ty can皿ot be overcome， Unless the round off errors are removed． For example dγ／dx−y−2θ”x， whose solutionγ＝θ一x with the initial valueアo−1 at x＝O tends to zero as x→oo， whie all other solutionsア＝．4ex十e−x，．4±O tend to ：ヒoo aS X→OO． We Can not pursue the solution y＝ e一ズby numerical integration by the round off ㎝ors． Some equations may have a movable singularity in a very small neighbourhood of x＝0． Consider dyfdx−一γ2， and the solutionア＝100／（1−100x）withアo＝100 at x＝0． If we eake the step size h＝1／10， then we are to obtain another solution． hthese cases（panial instabmaty）we must choose the step siZe， depending on the initial value of the equation， so that the num凱ical hltegration can．be properly done． Further the solution of some differential equation may increase or decrease sud− denly and rapidly i皿 the course of㎞tegration， and this can not be knoWn unless the solutjon reaches at the poi皿t． dうy ＋K（y2−1） 十y＝0，with large For example consider van der Pol’s equation dx2 dx K（say K−一 100）． The solution increases or decrease suddenty and periodically． In the above mentioned cases we can not solve satisfactorily the equation by numerical integration． In such cases we fail to obtain the requiエed solution， or we continue to obtam the solution beyond the singular point． Besides these theoretical d幻田culties round off errors can not be avoided， which make the numerical solution meaningless． Many of dec㎞al digits may notもe true even if it is wright on the theoretical ground． 13．Mo〔rre，s interval arithmetic method． One of the important methods‘‘Interval arithmetic method”is given by R． E． Moore．40｝By his Inethod some of the di伍一 culties mentioned above are removed， and the interval solution， ca11ed by him， which contahls the true solution， consider桓g round off errors， can be obta口led． But his method is not perfectly complete． He gives two methods on the same stand point． The equation studied by him bei皿g the、 equation of the fbmめノ／dx＝∫（γ）， he gives the region［γ1，γ2］， that is， A such thatγo∈∠1，アo being the血itial value．41｝ His丘rst method does not give any condition on．4， and choose h so thatγo十 F（A）⊂A． This shows that the solution ofめ7／dx＝∫（ア）withγ＝Yo at x＝O exists ・u・d・・n・≦・≦h・…nh・9iv・・Y（h）一・・＋TF（・・）・w・er…i・ap・・…A…er− mined by p． He does not give any condition onク． Hence this method is n㏄es− sadly not so precise． 40）R．E． Moore：Interval Analysis， Prentice−HaU， Inc．1966． 41）Region A may be the de6nition doma血ofノてア）or may be taken suMciently small． We use Moore’s notations，［］， F， A， Y（［0， h］）， etc．56
T．SHIMIZU
伍ss㏄・ndm・血・｛・・s・hb・h−
iε1γ。IK！lF（K）（γ。）1）1∫ 1）・・（ε1腎∼餅！）11K・ and givesノ鍾（o）such．that A（旬＝Yo十F（o）（γo）［0，乃］． Tben he tests whether γ（［0， h］）⊂ノ1（o｝or not． If．｝て［0， h］）⊂A（o｝， then γ（x）is the required solution． Otherwise putting A（1）＝］r（［0， h］）he tests whether｝て［0， h］，∠（1｝） ⊂A（1｝． 盾秩@llot． If Y（［0， h］，∠4（1｝）⊂∠4（1）then｝てx，∠4（1））is the required one． Otherwise puttillg．4（3L］r（［0，乃／2］，．4（1））he tCsts whether．4（3）⊂A（1｝or not． If．4（3）⊂A（1｝then ア（x，ノ｛（1））is the required one． Otherwise he contillues the sim且ar process as fbllows． ・4（1」γ（［0，乃］，・4（o））compared with A（o）； 、4｛3）＝｝て［0，乃／2］，∠4に｝）with A（1）； ノd（5｝＝ 】！（［0， 乃／4］， ノ4（3）） with ／1（3）； A（7Lγ（［0，乃／8］， A（5））with A（5）； A（9）＝Y（［0，乃／16］，．4（7））with zl（7， ……… A（2）一｝て［0，h］，メ（1｝）Wlth A（1｝； ．4（4）＝］r（［0，乃／2］，A（3｝）with∠4（3）； ．4（6」］r（［0，乃／4］，、4（5｝）with．4（5）； ．4（8L］r（［0；．乃／8］， A（7））with．4（7）； etc． The second method of Moore has some ambiguity． First： Let the equation be dγ／dx＝eア， then thC true solution withγ＝γo at x＝0 …一・・＋1・g1−≡，。・・nd・一・一・・i・am・va・1・・in・・1…y・ in倣・ca・e・一寄 fbr the method IK＝1． Ifεis taken as the least unit of the computer we take yo so thatεγo＞1． For such aアo， h is larger than x＝e−，・， which is血e si皿gular point of the true solution． Moore’s second method is tO take h so that the relatiye error of the numerical solution may b㏄ome sma11， but the convergen㏄of the Tayloピs expansion of夕ω with resp㏄t to乃， that is， the singular po垣t of the solution is not considered． Second：Moore assumes that the sequence A（1）， A（3｝，…， A（2n＋1｝terminates in a finite number of steps， that is A（2n＋1）⊂A（2カー1）fbr a 6nite〃．． Now consider the solutionγ＝1／（i−9ex）of the equationめノ／dx＝cylo withアo＝1 at x・＝O． We shall consider the method K＝1（the same result holding fbr any’1Gth method）・By乃＝Mεyo／10c『γ019＝∨願we consider C＝＝ 2，ε＝1／10， h＝1／20・Now ∠（°｝＝1＋［0，h］c，A（1）＝1＋［0，h］［0，c（1＋乃c）2］． Consider the upper bounds ao，01，α2，…of．4（o），．4（1），．4（2｝，…resp㏄tively． Then ao＝1．1，α1＝1．47373， a3 ＝＝1．50178，α5＝2．45884，α7＝101．9758， ag＝＝7．6006×1017◆ h・・n・・al…一・−1＋1。S．．、（…一・）・・， a・・＋・−1＋1。i、。（・…・）…㎞・・d・…h・・ a2n＋1＞α2ヵ＿1 it is su伍cient that（a2”＿1）9＞10・2カ．T・・h・w・hi・w・・㏄（・…）・〉・…〉｛織〉（1謬1°＞1・・・…if・・＞1・・…
On the other hand（α亙）9＞ax＞10・2K holds fbr K＝9． Hence the above inequality holds f（）r a皿〃≧9． Thus a2n＋1 monotonously and inde丘nitely垣creases and there is noπsuch thatα2潟＿1＞α2π＋1． The same result holds fbr the equation⑳μx＝θ， discussed befbre．・・…b・・h・皿…1・・1・・a・ ・…he eq…i・n，・nd・」要，w…eV…≧1．
Now， putting b＝V砺， ao＝γo十b， a1＝γo十bexp（b）〉γo十2．71，a3＝γo十（ゐ／2）exp（b exp（●）） 〉γo十7．57， α5＝γo十（b／4）exp（（b／2）exp（b exp（ゐ）））〉γo十488．1． In general a2n＋1＝YoCONT［RIBUTION TO THE THEORY OF NUMERICAL INTEGRATION
57
十（b／2りexp（a2銘一1）． By（6／2n）α2ヵ＿1＞a2”＋1， which is deduced丘om a2n−1十iog e＞〃log 2 −lo96十109α2カ＋1． Henceα2ヵ＋1＜a2ヵ一1 can llot o㏄urs fbr any n42） ． Though the五rst method of Moore is rigorous but it can not avoid the partial instability．…th・θq・・・・…⑳／血一・・w・・h・−i詰…一・・Ta・・［〒；，・］・・A，・n・
yo十xF（A）⊂．4， x＝1／3 andρ＝1／3． Now−1／100十（1／9）［−29／900，391／900］2＝［−0．0115，． 0．0109…］where 5＝［−29／900，391／90q］， which conta口1s a solution with positive initial values． Hence we see that the region／1 must be take皿so small that it does not． contain the posi葡ve value as the upPer end of A． Without choosing．4 and、ρproperly there may o㏄ur Fox’s partial hlstab伍ty． Let： ．4＝［−5，5］and the equation dγ／dx＝−Ly withγ（0）＝1． Th㎝乃is detemined by’ 1十5hL≦5，1−5んL≧−5，丘om which乃＝（4／5）」L． Choosmgク＝2 we have∫＝γo十［0，乃∫ρ］F（A）＝［1−2，1十2］＝［−1，3］，γ1＝γσ十・ h 万F⑲＝［1−6／5・1＋2／5コ＝［−1／5・7／51 ・ The true solution lying between O and l the upper solution and the lower one gives not pr㏄ise result at the next st6p and the following． On the other hand h＝（4／5）L． is so sma11， ifLis large． 14．On a n皿merical method． We sha皿give here a num｛）rical method fbUowing： Moore’s method，43）which can be used to integrate the equation without㎞owing・ the range of eXistence or the behavior of the solution． This corresponds to a、 mgthod． which chooses p orεof Moore’s so that it depends on∫（x，γ）and the initial value． ㎞the f（）110wing， in order to avoid rounding errors we must adopt rounding up’ computation fbr a皿calculation obtaining the upper numerical solution and ro皿d− ing down computation for all calculation obtai晦the lower numerical solution as Moore， s interval arithnetic method． Let the d田Terential equation be ・（51） 塞一f（・，・）
where∫（x，γ）is analytic with resp㏄t to x andアin some domain． 42） The proof that if A（2カ＋1）⊂A（2n−1）， then yo十xF（A（2”−1））contai皿s the true solution fbr へ 0≦ざ≦h／2n is not given by Moore． Let乃be the upper Hmit of乃such that yo十 ［0，h／2＃−1］F（A（2外一1））⊂A（2n−1）． Then if．4（2n＋1｝⊂ノ1（2n−1｝we have h／2＜万， this shows． that．4（2n＋1）contams the true solution， since γo十［0， h／2外一1］断てA（2カー1｝）contains the tlue one． 43）Moore assumes thatノてy）is to be a rational function ofγ． To calculate FてA）or F（S1）is tO calculate the maximum or the mini皿um of F血1．40r inぷ1． By interval arithmetic method we can obtaill the maximum or the minimum somewhat easier thall by or（㎞ary calculus， but the result is rough． If we can obttin max l f I or miniげl by ordinary calculus it is better than that obtailled 1）y血terval arithmetic method．：58 T．SHIMIZU We consider the solutionア（x）withア＝アo at x＝・O， and the numerical solution yヵ ：f（）rx＝＝xn， n＝0， 1，2，… and the step−size乃＞0． If f（x，ア）is de五ned in some domail1σso that the solution with some initial value ?wists locally and is unique， then we can integrate numerically the equation inσ． Though it is not known whether the solution exists throughout inσor not， we ぐcan do so at least in some interval O≦x＜η（也e solution going out of the domain σif the solution does not exist throughout）． In practiceア（x， y）is given by fbmula， fbr example，めノ／dx−・Vi＝lyFi’， and we must determine the domain in which the uniqueness of the solutioll is satisfied． For the above equation the linesア2＝1 are s垣gular lines where the uniqueness of Lthe solution is not satis丘ed． Hence when we血tegrate nume五ca皿y the above equa− ’・tion， we must determine the domain in which the solution exists and is unique， and we can not consider the solution outside the domah1レ1≦η＜1，0≦x＜oo． For the equation dγ／dx＝（x−1）／γfbr which the 1血eア＝O is the line where the エight hand side of the equation is not continuous， we must exclude the narrow band 」ア1＜η，befbre we bigin to integrate numerically the equation． Thus we must draw the curves on which the integration losses its meaning， when We bigin tO integrate． In the fbllowing we call such curves‘‘true boundaries．” ． The same integral curve is in different ways calculated by numerical integration， when it is expressed by different differential equations． For example the equation ．d2γ／dx2十ア＝o orめノ／dx＝z， dz／dx＝一γgives a sine−curve， but the equation（めノ／dx）2 ＝1一ア20rめノ燐＝∨T＝アgives only the part iyi＜10f the sine−curve． When y apProaches to 1， the uniqueness of the solution is broken， and numerical integration ．must be stopped．whereアreaches 1． S㎞ce the theory of difEerential equations does not give any conditions under which partial or inherent instability may apper， we can not show some methods by whi（）h sudh an i皿stabiiity may be avoided． But one of the reason why there may appear some di伍culties in the numerical integration is that the domain contains the region ．within which the sign ofめノ／dx or d2γ／dx2 changes． Hence we must divide the domain into such sub−domails， in which the sign of dy／dx and d2y／dx2 is constant． It is diMcUlt to divide a domain into such sub−domains in general， but in prac− ’tice it is not so dificult to do so． Now consider the sub−do皿ainsγ’＞0， y”＞0；γ’＞0，ア”＜0；γ’＜0， y”＞0；〆＜0， y”＜Orespectively， where〆≡dγ／dx，〆L d2y／dx2． In the fbllowing we call such curves which diヤide the domain into sub−domains 4‘モ窒奄狽奄モ≠戟@boundaries．，， Let us consider the五rst order differential equation dy／dx ・−fd〔x， y）． The critical lboundaries c6rrespond to the curves on which〆＝O or y”一一 O． III general we can not draw such critical bounda丘es in a required precision except
CONTRIBUTION TO THE THEORY OF NUMERICAL INTEGRATION
59
straight lines or simple curves． But we can remark such boundaries in the course・of numerical血te醐on．
Consider the domain D， containing the initial point（0， yo）， bounded by true bound− aries． Ih that domain we draw all critical boundaries． Further consider the band−domains B1，β2，…which are bounded by two curves． ＞parallel to the critical boundaries on both sides of width 2ηor bounded by curves・ paranel to the true boundaries of widthη（ηbehlg sma11）． We ca11‘‘ihn〔）r domain”the domain excepting all the band−domains of two kinds． above． ． （αッo） reglon critical ヅ＝0 Fig．1 critical boundary ary of the band−domains containing the critical boundary or of the band−domains． bounded by the true boundaries‘‘the side−curve of the second kind．”・…xamp1・・D・a
秩Eurv・・」」∫（…）一・・nd・〃一纂＋●・（・h・・i・・
critical boundaries）． Let ab and cd be two of them．The above mentioned band−domam，
D＊＊，side curves etc． are shown㎞Fig．2． We shall obtain in the f（）nowi皿9伽o numerical solutions， the upper solution and the lower one， so that the true solu− tion of the equation lies between them． （x，γ， When the numerical solution reaches the side−curve of 1）1＊＊， two cases o㏄ur． If the upper solution or the lower one m㏄ts at（h1， Yi）the side curve of the （x，y。， 丘rst㎞d， then we consider the domainD2（h、≦x≦X，，ア、≦γ≦壬）and a part 1）、＊＊ ・
of D2 belonging to 1）＊and is connected． Fig．2 Thus we continue the same process as done f（）r the part D1＊＊． If the upper solution or the Iower one meets at（h1， Yl）the side curve of the s㏄一 We call“the趾st region”apart D＊of’ the imler domain which contains（0，γo）and is connected． Now consider the domai皿1）1（0≦x≦X，， γo≦ア≦γ）and a part 1）＊＊of D1， bdonging・ to 1）＊and is connected． We call a part of the boundary curves of’ D＊＊which belong to the boundary of 1）1 ‘‘ 狽??@side−curve of the 丘rst kind” and a． part of the boundary curves of D＊＊which． belongs to the boundary of the band−bound−・・ ． αB1
6 side curve of iαY） tbe lst kind C D 1）＊＊B2
7 （O．yo） 1 ﾜ｛e2認「1…。3f∂ ’1
60
T．SHIMIZU
ond kind， then we do not consider 1）2， and we consider a domain D＊＊＊which is bounded by the side curveγ， a curve parallel toγand two straight lines x・＝h1， x＝X，parallel to theアーaxis． ’ Thus we continue the same process as done fbr the part D1＊＊． Since the upper solution differs丘om the lower one by a small quantityδ， bqth the upper one and the lower have their end points on the side curve of the same ．㎞d． We cail‘‘the second region”the domain which is a part of the口mer domain and is conn㏄ted and lies op the other side of the critical curve （to the first region） and i皿to w］bich the solution enters， traversi皿g the critical boundary． If the numerical solution meets the side curve of D2＊＊at（h2， Y2）， and the side− curve belongs to the first㎞1d， then we consider D3（h2≦x≦X，，ア2≦y≦Y，）and 1）3＊＊of 1）2 belonging to D＊and is connected． We continue the same process as I）efbre and so on． If the side−curve belongs to the s㏄ond kind or if the numerical solution meets 裳 the side−curveγ20f D＊＊＊， we consider a band−domain 1）z＊＊＊of width smaller than the fbrmer， which is bounded byγ2， a curve parallel toγ2 and two straight Hnes parallel to theγ・axis． We continue the same process as befbre． In this case the band−domains D3＊＊＊，1）4＊＊＊，…beρome narrower and narrower， and the step size becomes smaller and smaller． If the number of d㏄imal places is fixed， we can not choose h so small， and the partial instabiHty can not be avoided． If we can increase the number of dec㎞al places， then by choosing smaU h， we can decide whether the upper or the lower solution remains in the丘rst region or traverses the critical boundary and enters into the second region． In the second region we cont垣ue the same process as done fbr the first region． For the system of equationsめノ／ゐζ＝∫（x，ア）（7 and f being vectors）we must consider a四rve in n−dimensional space． It is not easy to draw the region above mentioned． For example we consider a system of two equations 書∫（・，・，・），農一・（・，・，・）． 」〆＝Oor z’＝O respectively， and F（x，γ，z）＝O or G（x，y，z）＝O on which y”＝O or z”＝・0． ア（x）and z（x）are the proj㏄tions of the space curve． When the solution intersects the above surfaces， thenア’（x）， z’（x）， y”（x）or z”（x） vanishes． If we can see that the solution approaches the above surfaces we can do similarly as fbr the first order equation． In some cases we can see it analytically and the instability can be avoided．CONTRIBUTION TO THE THEORY OF NUMERICAL INTEGRATION
61 ‘ For example：For the equationsめノ／dx＝z， dz／めζ＝一アif yωor z（x） approaches to the z＝O plane or to the y＝O plane， then y（x）takes the maximum value or the point of inflexio皿． We can not draw c亘tical boundaries fbr the system of equations， and even if the equation is of the丘rst order， we can not draw critical curves pre（iiSely． Hence in practice， when we are to integrate the equation or the equations we do not draw the regions above mentioned， and we shall proceed as follows． In the course of numerical integration we calculate on the other hand the values off（x，y，z），9（x，y，z），F（x，ア，z）andG（x，γ，z）． If one of them at least becomes small in absolute value， then we must diminish 『h，and integrate the solution Carefully as mentioned above． By this process we can follow approximately the true solution， when the numerical solution goes out of the first region and into the s㏄ond one． ’ 15．Ghoice of the step size． First we must determine the step size of numerical int・9r・ti…fth・diffe・enti・1・q・・ti・n i…d・・t・qbゆth・upPer s・luti・n・nd th・ lower one． 、 ． It has never b㏄n discussed that how the step size must be determined in order to avoid instability in the course of llumerical integration． In the books fbr numerical analysis the fbllowhlg remarks w巡be fbund． For one step method consider the difference of the value Obtained with the step s元ze h and that obtained 1）y two steps with the step size乃／2． For predictor−corrector method consider the difeerence of the value obta㎞ed by the predictor and that obtained by the￠orr㏄tor． If the difference becomes larger than a fixed value， then it is recommended that the ．step size may be decreased． in order to obta口1 precise result． ． Another method is given by Moore as the second method of the interval method．・…t・p・セ・h−
iε1夕。IK！1ア（x）（o）1）11K ・…h・K−・h・・d・・m・・h・d・…9i・・…h㎞・ For the ditferential．equationめノμx＝fd〔y）consider the Taylor’s expansion y（h）＝ γo十hyo’十… 十1／π！乃『γo（π）十… at x＝0， the radius of convergence ア bemg given by ≡（n！Yo｛n））1㌦・・一・…（・）・ For elementary function f（ア）of Y， fbr examples，ノ」γ2，θア，…we have yo㈲÷ 刀！γ♂＋1，÷（n−1）！enア・，…． H・・ce・h÷（9−’i−hII：，1’：1）11外・÷（ε1窪）1ノ鰺・…・ ’lherefore ifε1γo l＜1， then h is not larger than the distance of the s垣gularity of the solution丘oln x＝0． ‘ ・噸rsu・h・n剛・・i・n d・／・・一．f（・）， w…e・・ω！・ぷ田・eS・t眺（y荒｝）1”62
T．SHIMIZU
＜（（lti2 iyt（pt＋6＞！）1畑＋1’＜…＜÷， and we can not take the s・ep sセe h〒（〒｝織｛…÷）11x f（・r a 丘xed。K since h may be larger than the distance of the shlgularity of the solution 丘omx＝＝O． We．shall give here some exact method to determine肩br the equation dy／dx− f（x，γ）． In order to choice a suitable乃we must remark the following thr㏄conditions： The五rst condition：We must not take such a step siZe， by which one integrates beyond the shlgular pomt of the true solution． The s㏄ond condition：There may appear no hlstab伍ty， if possible． The third condition：Numerical solution has to approximate the true solution in high precision， if possible． First condition：乃can not be larger than the distance of the（movable or fixed） singular point of the solution加m x＝x‘． To obtain such a h there are two methods， one of which is given by the radius of convergence of Taylor’s expansion of the true solution at a point x＝xi， and the other of which is given by the maximum and the mhimum of l∫（x，γ）l in some domain considered． We must have by Tay10r’s expansion h＜司証（n！Yi｛n））11n whereγ‘ωis the value ofγ（勃）（x）at x＝xi． On the other hand let M＝max∫（x，ア）， m＝mini∫（x，γ）where U（0≦x≦X，アo≦ x，ノ∈σ 潔，7∈σ y≦］r）ifγo’＞O and （ノ（0≦x≦X，夕o≧ア≧］r）ifγo’＜0． Then Yo十hM≦Y if γoノ≧0 アo十乃〃1≧】ノ if γ’≦0． For the system of equations dy／ゐc rf（x，ア）， whereγand∫（x，γ）are v㏄tors， we must take the least of a且the rad五〇f convergence of the component fun（ン tions Oノ： On the other hand we consider the ’ maxima iM and the minima湖of the com− ponents of the v㏄tor functio，n and letハイ＝max‘ハ4， m＝mini湖． Thus we have γo十hM≦｝二where Y is the minimum of iY if iアo’＞0 γo十乃〃≧Y， where Y is the maXimum of iγ if iアo’＜0． ・ The radius of convergence being㎜iquely dete㎜ined if the solution with the血一 tial value is fixed， the maxima and minirha give different values of h when i Y are chosen properly．The second condition：We must determine h so that there may not o㏄ur
instability．44） 44） S血ce mherent i lstability does not o㏄ur fbr the stable state of the solution in prac− tice，鉦is not so important血practical analysis． ’CONTRIBUTI（）N TO THE THEORY OF NUMERICAL INTEGRATION
63
For the first order equation we choose乃so thatγo十乃M is contained m『the inner domain in§14． Ifγ。’＞O then we consider the domain［1（0≦x≦X， y。≦y≦γ）in which Max∫−Mini∫is sma11． By this process we can obtain h so that no partia1 ロ instab狙ty may o㏄ur at，the丘rst step丘om x＝0． Next we choose h so that the upper solution and the lower one do not enter into aband−domain bounded by critiCal boundaries h1§14． If we can choose h as small as we please（that is， we can increase the number of dec㎞al pla㏄s mdefinitely）then it wi皿be sure that partial instabdity and others Wi皿 be avoided． l In practice it is not so， and hence when the upper solution or the lower one enters into a band−domain after n−steps it may happen that the instab丑ity o㏄urs． Thus the partial instability can not be avoided fbr the computation． In order to decide whether the instabm ity has oocured or not， we proceed to obtain the upper solutions and the lower ones fbr many steps． If the d醗r㎝ce of the upper solution and the lower one has b㏄ome larger， then we see that the instability has oo｛ ured． Inherent mstability and partial one can not be distinguished for’numer− ical analys姪． The same phenomena occur fbr the solution whiCh varies its direction rapidly or which osc狙ates hl high−frequency． The fbmler case is， fbr example， the solution of van der Po1’s equation〆ノ十K（γ2−1）〆十γ＝O fbr larger K（say K＝100）． The latter case is， fbr exanlple， y＝sh11／（1−x） as x→1． If the difference of the upper solution and the lower one becomes large， then we take a smaller乃， as the first step from x＝0， than that was used befbre， and then we proc㏄ed as befbre selpcting乃1；h2，…． Of course it may happen that．the sequenCe of the step sizes h1， h2，…b㏄omes less than the least unit of the computor． Then the process must stop． if the number of steps increases then the round up or down errors also increases， but it will give a better result if the皿umber of decimal places can be increased． For the system of equatiolls we can not draw critical bo皿daries， or band−domains， but if we can see the num頃cal solution approaches the above critical surfaces， the ’ partial instabiity w皿be avoided． It will be done in some cases， because we are calculating the values of the crit− ical sur白㏄s in the course of numerial垣tegration of the equation as mentioned in the last part of§14． Thus we can cont血ue the numerical integration unt皿the numerical solution approaches the true boundary． The third condition：We must determine h so that the numerical sOlution ’may be obtahled vely near the true one． We shan here adopt Taylor’s expansions method of order K≧1 （K＝ 1 being EUler，s method）御
T．SHIM工乙U
1．∫、・αサ・（・）吻’（・）＋芸・〃（・）＋…＋姜；ノ…（・）・ Now’we suppose the order of aocuracy may be 10−”1． ＿ 、 ・ir・t we ch…e・…u・h…g・h・・芸・…（・）−1・一・… Next consider the sign ofア’（0）andン”（0）． Ifア’（0）＞0， y”（0）＞0（other、 cases「wm be s㎞皿arly done）． ”・，Letσbe a doma㎞．in whk海、the sign ofγ’ωandγ”（x）is constτant． ∵Forヰhe first Order eqqation it has been disCussed． in．§14・ For the system of m equations consider the domain（朋十1 dimensional domahl） 0≦x≦x，iy。≦輌y≦‘γi if‘γ。’＞0；．げ。’−0，げ・”＞0； ．Ui・≦・≦X，、。。≧i7≧、Y、 if、。。’＜・、、。。⊆・，、・。〃＜・・ ．where《γ，ゴγ，…are the components of．the v㏄torγ（the sign ofεγ’， iy’・ξγ”・．JY” be㎞g collstant虹1こ九）． 、 」 Let ， 、 ∵．D M・−Max∫（x，y）・m・−Mini∫（x・γ）i・U・・
Now we choose
〃o＋‘乃1’Ml＝iYl fbr each component． ‘ 」γ。＋」乃2’殉r】三 Let h1，l be the sma皿est of乃o．乃1’， h2’；乃1’，乃2’behlg． the mhlima of《乃1’，￡乃2’ Iespeetively． 一． ・ ． ’ （Since there誌nb ambiguity we shall omit the prefix show㎞g the components of the vector） 、 Then we llave fbr each prqlection of、 the solution on the compone皿t−Plane アo十x〃11≦ア（x）＝3アo十零ゾ（θx，．プ（θx））≦yo十xハ4㌦ fbr O≦x≦乃1，1， wh・・eγωi・th・tm・・61uti・n a・d th・・ight h・nd・id・・f． y（x） iS・b楓in・d b測ea・ valUe theor㎝． Now we apPly．Taylor’s expa皿sio皿method to obtam也C nu亘e工ical solutやn 9￡γ〈x）・ Thus we hatie the prqlection of two solutions（也e upper solution and the lower ＜）ne）on the component−pla且e， ・・＋乃・’＋…＋（xk−1k−1），・…一・・＋li妬ω ・ 一・’＋…＋（差∴，・…一・・＋芸m・働where Mi㈲，〃1（k｝denote the’． maximum and minim面ofア㈲mU： respectively．
The proj㏄tion of the． tme solution lying betw㏄n the projectio皿of the two solutions on the component−Plane． Letγ1＝maXimum of（the projecti（）n of the upper solution−that of the lower solu− i＝1， ゾ em tion）at X＝乃1，1． Ifγ1≦10一殉， then we take乃1，1 as the first step sセe． ・ 一 「］f．γ三＞10一餌1−then W7e take the domam（〃2十1 dimensional domah1） ・． U，（0≦x≦X2， iYo≦‘ア≦i】「i）・r（0≦x≦X2，」γ。≧∫ア≧iY，）， X2＜−、， Y2＜Y・ ot Y2＞γiCON’1 RIBUTION TO THE THEORY OF NUMERICAL INTEGRATION
6．5 alld choose乃1”， h，”the m㎞㎞a of ihl”，《乃2”so that ・ 一・ ぎア・＋乃・”M2−・Y・． ，ア。＋，h、”m、r｝ち’』h−’− where M2， M2 the ma∼dmum， the minimum． of∫垣σ』respectively． We take hi，2 the smanest ofゐo， h，”； h2”．二 ・二Considering TqylorヴS expansfon method −− tt・・櫛臼…＋（xk−1k−1），・…一・・＋＞iM・㈲輌・’＋…＋蒜，・…一・・＋li垣・ω 1
’where妬㈲， M2㈲the maxim㎜， the ．minimum of y（k） in［12 respectively． Letγ2＝the maximum of（the pr（）jectiOn of the upper solution−that of∫heエlowerS？耳uti・n）・t x一乃…三． ． r
Ifγ2≦1ぴ餌1， then we take乃1，2 as the丘rst step sセe． ・・Ifγ・＞10 ・， th・n we cC・ti・…画・・斑・㏄S・．a・b・fo・e． The domainsこ1，，こ13，…becomin9．smdlldr and the max㎞a（minima）ルt2，妬，… イ〃12，〃23，・、・・）d㏄reasing（increasing），γ2，γ3，… decreases indd〔initely． r】「Cthe number of d㏄hnal places is larger than〃11， we haveγ♪so thatγρ≦10一傷1 hnd we take乃1，♪as the first step size．（lf．K is su伍cielltly large， then乃o may be su伍ciently small，もut if、K＝1，0r’2thell ho is not so sma皿and乃o，乃1，、h2 are of ・compatible size）．子 1． ． Let the prqj㏄tion of the two solutions45）（ttie upper solution and the lower one） 乃1，ρ頁 γ。＋乃、，〃。’＋…＋ 妬（KL五 K！ ・・＋・・，…’＋…＋讐殉…一控 Considerア’（乃1，ρ）＝max∫（h，，♪，γ）fbr些≦γ≦五，どノ（h1，♪）＝mini∫（乃1，♪， y）負）r∼2≦γ ≦完． A㏄ordhlg as the pr（）jection of冗ノ＞O or that ofど2’＞0， we consider ぺ U・（h1．ρ≦X≦X，ど1≦ア≦】rl），些≦｝コ・rU，（h、，♪≦X≦X，〕㌃≧y≧竺），五≧五． lSin㏄寅一ど1≦10一協1，五is of the same sign withど1 in general but if the sign of五 孤dど！is different we must take the order of aceuracy 10一郁smaller than the given ．10−”1． ；Let Mi＝maxf（X， y），〃1、−minif（X，γ）in U、． ・Choose h1’， h2’the minima of《乃1’，‘乃2’resp㏄tively so that ㌶＋ih・’M、一菰，辺＋ih2’〃、 r五． C・…eh・一㎡・・（h・…’・h・’）w・er…一（，、・緩，，））1’κ・ Then the proj㏄tion of the solution ム＋X〃1・≦yω＝y（h、，♪）＋xf（θX，夕（θX））≦五＋Xハ4 fbr h、，ρ≦X≦h2． 45） In the s㏄ond step，茄， Yi，〔九，ハ41， Ml， hl’， h2’etc． are different from those in the first step， but are expressgd with the same letters．66
