Advances in Mathematical
Sciences and Applications
Vol. 26, No. 1 (2017), pp. 29-38

GAKKOTOSHOTOKYO JAPANNUMERICAL LIMIT AND ITS APPLICATION
TO A BLOW-UP PROBLEM RELATED TO DEFAULT RISK
Yuya Sasaki
Komi High School,
1006-2 Chiyosato, Komi, Minamisaku, Nagano, 384-1105, Japan(sasaki-y@nagano-c.ed.jp)
Hiroko Soutome
Japan Woman’s College of Physical Education,
8-19-1 Kitakarasuyama, Setagaya-ku, Tokyo, 157-8565, Japan(sotome.hiroko@jwcpe.ac.jp)
Hitoshi Imai
Faculty of Science and Engineering, Doshisha University,
1-3 Tatara Miyakodani, Kyotanabe, Kyoto, 610-0394, Japan(himai@mail.doshisha.ac.jp)
andNaoyuki Ishimura
Faculty of Commerce, Chuo University,
Hachioji, Tokyo, 192-0393, Japan(naoyuki@tamacc.chuo-u.ac.jp)

Abstract. In the paper, we propose the method of numerical limits, which consist of theextrapolation and the bounding transformation. We apply these numerical limits to thecomputation of the blow-up time for a model ODE system of interacting firms. Numerical
implementation shows satisfactory results.

————————————————————
Communicated by Editors; Received January 8, 2017
This work was supported by JSPS KAKENHI Grant Numbers 16K13774 and 16H02155.
AMS Subject Classification: 65-05, 65L05.
Keywords: blow-up, numerical limit, infinity, differential equation.

29

30

1

Introduction

In the paper, we are concerned with a numerical treatment of the limiting behavior of
blowing-up solutions. In particular, we apply our method to the following system ofordinary differential equations (ODEs), which is introduced in Ishimura and Nakamura[5][6] as a simple model dynamics of comercial firms:
0
u (t) = − α u(t) + α v(t) w(t)v 0 (t) = − β v(t) + β w(t) u(t) 0w (t) = − γ w(t) + γ u(t) v(t)(1)

where αβγ are positive constants.
As is well known, blow-up phenomena for differential equations have been investigatedextensively so far. We recall the seminal work of Fujita [2] for the semilinear heat equation.
Here we mean that a function u(t) = u(t, ·) blows-up if lim ku(t)k = ∞ with k · k beingt↑Ta suitable norm. T (< ∞) is called a blow-up time of a function u(t). It is shown thatthe solution to (1) may blow-up or breakdown within finite time [6], that is, according toinitial conditions there corresponds 0 < T < ∞ such that |u(t)| + |v(t)| + |w(t)| → ∞ ast ↑ T . In [6], the blow-up time is interpreted as a default time of firms.
To illustrate the novelty of our nethod, we consider the following two cases of (1) asshown in Table 1.
Table 1. Typical two cases of (1).
(α, β, γ)(u(0), v(0), w(0))

Case 1( 10, 10, 10 )( 10, 10, 10 )Case 2( 1, 10−6 , 100 )( 10−6 , 100, 1 )It is noted that Case 2 treats extremely different values of the coefficients and the initial
data.
Numerical solutions of (1) are depicted in Fig.1, where the Euler method with the1variable temporal mesh size [7] ∆t = min, ∆t0
is used. Here ∆t0
max{|u|, |v|, |w|}
denotes some fixed constant specified in each implementation. Numerical computation is
performed in double precision. Obviously, in Case 1 u(t) = v(t) = w(t) holds theoreticallyand numerically; indeed, three profiles of u, v and w in Fig. 1 (1) (Case 1) are seen to
ovrelap each other. Fig. 1 (1) (Case 1) suggests that the blow-up times for u, v, w are
the same. On the other hand, Fig.1 (2) (Case 2) indicates that the blow-up times may
be different, while this is not true according to the theorem in [6]; for this system, it isproved that blow-up times for u, v, w are the same.

31
10000

10000

u
v
w

8000

6000

6000

4000

4000

2000

2000

0
0

u
v
w

8000

0.002

0.004

0.006

0.008

0.01

0.012

0
0

0.002

0.004

0.006

t

0.008

0.01

0.012

0.014

0.016

t

(1) Case 1.

(2) Case 2.

Fig.1. Numerical solution profiles (∆t0 < 10−5 ).
Our numerical method gives the numerical blow-up time Tnu as shown in Table 2; here
Tnu is regarded as the blow-up time when ∆t < 10−16 .
Table 2. Numerical blow-up time Tnu .
∆t0Case 1Case 210−50.0105434262 · · ·0.0159578799 · · ·10−6

0.0105370056 · · ·0.0159632635 · · ·10−7

0.0105361214 · · ·0.0159638484 · · ·We should remark that this numerical method does not apparently guarantee that blowup times for u, v, w are the same. On the other hand, our method yields the same blow-up
time even numerically.
The numerical computation of blow-up problems is challenging, since the blow-up
phenomena involve limit and infinity. To list up some literatures, we recall the workof Nakagawa [7], where an adaptive temporal mesh control as shown in the above is
proposed to compute the blow-up time, and it is shown that the numerical blow-up time
converges to the exact analytical one. Recently, Cho [1] proposes a nice scheme with theuniform temporal mesh size and shows the convergence property as well. In his scheme
the difficulty of infinity is removed theoretically, and the difficulty of limit is partiallyremoved by taking the temporal mesh size be smaller.
In the present article, we propose a numerical method which enables us to remove thedifficulty both of limit and inifinity. The method consists of the extrapolation and thebounding transform [4]. It is applied to (1). Our interest is not only in the numericalcomputation of the blow-up time, but also in the possibility that numerical results may
show that blow-up times for u, v, w are the same.

32

2

Our methods

To remove the difficulty that the valueof a function u becomes infinite, we consider the
next bounding transformation [4] of u into ũ.
!
√ũ− 1 + 1 + 4u2
u=
.
(2)
ũ =
2u
1 − ũ2
By this bounding transform, there follows the onto-mapping between the intervals of u
and ũ.
ũ ∈ [0, 1] ⇐⇒ u ∈ [0, ∞].

(3)

There are two ways for the application of the bounding transform in numerical computation. The first way is to numerically compute the transformed equations, which is
known to become more complicated than the original equations. The second way is toperform the bounding transform to a numerical data, which are obtained from original
equations. The second way is simple enough to use it here. We should remark that in
the current numerical computation, a standard double precision does not work effectively.
For example, due to the loss of trailing digits, ũ = 1 corresponds to u ' 1016 in double
precision. Therefore. multiple-precision arithmetic is necessary.
To remove the difficulty of limit we use the extrapolation. For example, we fix the
temporal mesh size ∆t and then the extrapolation is employed for numerical computation:
lim . The details of our procedure will be explained in the next section.
∆t→0

3

Numerical results

The forth order Runge-Kutta method is used for the discretization. The temporal meshsize ∆t is fixed. Numerical computation is performed in multiple precision (100 digits)
by using exflib [3] which is a compact and fast library.
Firstly, we show numerical results for Case 1.
10000

1.00

u
v
w

8000

ũ
ṽ
w̃

0.99

6000

0.98

4000

0.97

2000

0.96

0

0.95

0

0.002

0.004

0.006

0.008

t

(1) Original profiles.

0.01

0.012

0

0.002

0.004

0.006

0.008

0.01

0.012

t

(2) Profiles of transformed solutions.

Fig.2. Numerical solution profiles in Case 1 (∆t = 10−5 ).

33
In this case u(t) = v(t) = w(t) holds theoretically and numerically. It is seen that threegraphs in each of Figs. 2−5 overlap each other. The original profiles of u, v and w
obtained by the forth order Runge-Kutta method with ∆t = 10−5 are shown in Fig.2 (1).
Of course, profies are almost same as in Fig.1 (1). Then we apply the bounding transformto these numerical data of u, v, w. Profiles of transformed solutions ũ, ṽ and w̃ areshown in Fig.2 (2). It is interesting that profiles of ũ, ṽ, w̃ are almost straight lines. To
get better straight lines, we consider a narrow interval [0.0100, 0.0105] which is close to
the blow-up time. Numerical data of ũ, ṽ, w̃ in this interval are plotted and regression
lines are exhibited in Fig. 3. The t−coordinate at the intersection of two lines, namely
the intersection of ũ = 1 with the regression line for ũ, means the numerical blow-up timeof u for ∆t = 10−5 . We denote this time as Tu (∆t = 10−5 ). We obtain Tv (∆t = 10−5 )and Tw (∆t = 10−5 ) similarly. We note once again that in Case 1 these three numericalblow-up times are the same.
1.00000.9995ũṽ
w̃

0.99900.99850.99800.99750.99700.01000.01010.01020.01030.01040.0105tFig. 3. Regression lines in Case 1 (∆t = 10−5 ,t ∈ [0.0100, 0.0105]).
In Fig. 4, Tu (∆t), Tv (∆t) and Tw (∆t) are plotted for various ∆t, and regression lines areshown. The vertical intersection with the regression line for ũ means lim Tu (∆t)(≡ Tu ).
∆t→0

From numerical data in the interval [0.0100, 0.0105] we obtain that Tu = Tv = Tw =0.01053650 · · · .
0.01053651

Tu (∆t)Tv (∆t)Tw (∆t)0.01053651

0.01053650

0.01053650

0.01053649

0.01053649

0.01053648
02x10 -6

4x10 -6

6x10 -6

8x10 -6

1x10 -5∆tFig. 4. Numerical blow-up times in Case 1 (t ∈ [0.0100, 0.0105]).

34
Similarly, from numerical data in the interval [0.01050, 0.01053] which is closer to theblow-up time, we obtain that Tu = Tv = Tw = 0.01053605 · · · (Fig. 5). Please compare
this with Tnu = 0.0105361214 · · · in numerical computation with variable ∆t (Table. 2).
0.01053605435

Tu (∆t)Tv (∆t)Tw (∆t)0.01053605435

0.01053605434

0.01053605434

0.01053605433

0.01053605433
05x10 -8

1x10 -7 1.5x10 -7 2x10 -7 2.5x10 -7 3x10 -7 3.5x10 -7 4x10 -7 4.5x10 -7 5x10 -7∆tFig. 5. Numerical blow-up times in Case 1 (t ∈ [0.01050, 0.01053]).
We remark that the exactness of our numerical results depends on the interval of t.
We apply the same proceedure to Case 2. The original profiles of u, v and w obtainedby the forth order Runge-Kutta method with ∆t = 10−5 are shown in Fig. 6 (1). Of
course, profies are almost same as in Fig. 1 (2). Then we apply the bounding transformto these numerical data of u, v, w. Profiles of transformed solutions ũ, ṽ and w̃ areshown in Fig. 6 (2).
1.0

10000

u
v
w

8000

ũ
ṽ
w̃

0.8

6000

0.6

4000

0.4

2000

0.2

0.0

0
0

0.002

0.004

0.006

0.008

0.01

0.012

t

(1) Original profiles.

0.014

0.016

0

0.002

0.004

0.006

0.008

0.01

0.012

0.014

0.016

t

(2) Profiles of transformed solutions.

Fig. 6. Numerical solution profiles in Case 2 (∆t = 10−5 ).

35
To get better straight lines we consider a narrow interval [0.0155, 0.0159] which is close
to the blow-up time. Numerical data of ũ, ṽ, w̃ in this interval are plotted and regression
lines are shown in Fig. 7.
1.000000

0.99980.9996ũṽ

0.999999

0.99940.99920.999999

0.99900.99880.999998

0.9986
0.99840.999998

0.99820.99800.999997
0.01550.01555

0.01560.01565

0.01570.01575

0.01580.015850.01590.99780.01550.01555

0.01560.01565

0.0157tt(1) ũ.

(2) ṽ.

0.01575

0.01580.015850.01591.00000000

w̃
0.99999995

0.99999990

0.99999985

0.99999980

0.99999975
0.01550.01555

0.01560.01565

0.01570.01575

0.01580.015850.0159t(3) w̃.

Fig. 7. Regression lines in Case 2 (∆t = 10−5 , t ∈ [0.0155, 0.0159]).

In Fig. 8, Tu (∆t), Tv (∆t) and Tw (∆t) are plotted for various ∆t, and regression lines
are exhibited. From numerical data in the interval [0.0155, 0.0159] we obtain that Tu =0.01596064 · · · , Tv = 0.01597657 · · · , Tw = 0.01595765 · · · . The first 3 digits are the
same.

36
0.0159766

0.01596066

Tu (∆t)0.01596066

0.0159766

0.01596066

0.0159766

0.01596065

0.0159765

0.01596065

0.0159765

0.01596065

0.0159765

0.01596065

0.0159765

0.01596065

0.0159765

0.01596064

0.0159764

0.01596064

0.0159764

Tv (∆t)0.0159764

0.01596064
02x10 -6

4x10 -6

6x10 -6

8x10 -6

1x10 -50

2x10 -6

4x10 -6

6x10 -6

8x10 -6

1x10 -5∆t∆t(1) Tu (∆t).

(2) Tv (∆t).

0.01595772

Tw (∆t)0.01595771

0.01595770

0.01595769

0.01595768

0.01595767

0.01595766

0.01595765
02x10 -6

4x10 -6

6x10 -6

8x10 -6

1x10 -5∆t(3) Tw (∆t).

Fig. 8. Numerical blow-up times in Case 2 (t ∈ [0.0155, 0.0159]).
Similarly, from numerical data in the interval [0.0159, 0.01596] which is closer to the blowuptime we obtain that Tu = 0.01596394 · · · , Tv = 0.01596397 · · · , Tw = 0.01596386 · · · (Fig.
9). The first 5 digits are the same. Please compare this with Tnu = 0.0159638484 · · · innumerical computation with variable ∆t(Table. 2). Numerical results suggest that solutions blow up at the same time.

37
0.0159639437

0.0159639796

0.0159639436

Tu (∆t)Tv (∆t)0.0159639795

0.0159639436
0.0159639794
0.0159639435
0.0159639793

0.0159639435
0.0159639434

0.0159639792

0.0159639434

0.0159639791

0.0159639433
0.0159639790
0.0159639433
0.0159639789

0.0159639432
0.0159639431
01x10 -7

2x10 -7

3x10 -7

4x10 -7

5x10 -7

6x10 -77x10 -78x10 -7

9x10 -7

1x10 -60.0159639788
0

1x10 -7

2x10 -7

3x10 -7∆t

4x10 -7

5x10 -7

6x10 -77x10 -78x10 -7

9x10 -7

1x10 -6∆t(1) Tu (∆t).

(2) Tv (∆t).

0.0159638664

Tw (∆t)0.0159638663

0.0159638662

0.0159638661

0.0159638660

0.0159638659

0.0159638658
01x10 -7

2x10 -7

3x10 -7

4x10 -7

5x10 -7

6x10 -77x10 -78x10 -7

9x10 -7

1x10 -6∆t(3) Tw (∆t).

Fig. 9. Numerical blow-up times in Case 2 (t ∈ [0.0159, 0.01596]).

4

Conclusion

In the paper, a blow-up problem governed by the system of ordinary differential equations
is analyzed numerically. The system is related to the default risk and various properties
of solutions are theoretically clarified. The interesting property from a view point of
numerical analysis is that solutions blow up at the same time. Traditional numericalmethods are not suitable for ascertaining this property. We have here, on the other
hand, proposed the method of numerical limits which consist of the extrapolation and thebounding transform. We apply this method to the problem of the numerical computationof blow-up phenomena. Numerical implementation is performed in multiple precision and
it is realized that our procedure works efficiently. Numerical results successfully indicatethat the blow-up time are the same.

38

References[1] C.-H. Cho, On the computation of the numerical blow-up time, Japan J. Indust.
Appl. Math., 30(2013), pp. 331-349.
[2] H. Fujita, On the blowing up of solutions of the Cauchy problem for ut = ∆u + u1+α ,J. Fac. Sci. Univ. Tokyo Sec. 1 A Math., 16(1966), pp. 105–113.
[3] H. Fujiwara and Y. Iso, Design of a multiple-precision arithmetic package for a 64bit computing environment and its application to numerical computation of ill-posed
problems, IPSJ J., 44(3) (2003), pp. 925–931.
[4] H. Imai, Some Methods for Removing Singularity and Infinity in Numerical Simulation, GAKUTO Internat. Ser. Math. Sci. Appl., 23(2005), pp. 103–118.
[5] M. Inoue, N. Ishimura and M. Nakamura, Numerical Study on the Systems of Nonlinear Ordinary Differential Equations for Default Risk Model, in IEEE 17th International Conference on Industrial Engineering and Engineering Management (IE&EM),(2010), pp. 1066–1069, availabe at http://ieeexplore.ieee.org/document/5646439/.
[6] N. Ishimura and M. Nakamura, Blow-up Behavior of a System of Nonlinear Ordinary
Differential Equations. J. Res. Inst. Sci. Tech., Nihon University, 2013 (2013), pp. 1–
8.
[7] T. Nakagawa, Blowing up of a finite difference solution to ut = uxx + u2 , Appl.
Math. Optim., 2 (1976), pp. 337–350.