Asymptotic Methods of ODEs: Exploring Singularities of the Second Kind

where A is an n n matrix with elements depending on t and w t is a column matrix whose elements are unknown functions. We will suppose that A has analytic coefficients (near infinity) and can be written as a series A j 0 t j A j (convergent for large t) for constant matrices A j where A 0 is non-trivial. (Such a singularity at finite t may be turned into a singularity at infinity by a change of variables.) This type of singularity for r 1 is called a singularity of the second kind (also known [4] as an ‘irregular’ singularity). Differential equations of the form t r j j o n a j t d n j d t n j y 0 a j ’ s bounded for large t>0 and a 0 1) can also be written in the above general matrix form by setting Asymptotic Methods.nb 8/5/11

The Mathematica Journal volume:issue © year Wolfram Media, Inc.

Differential equations of the form t r j j o n a j t d n j d t n j y 0 a j ’ s bounded for large t>0 and a 0 1) can also be written in the above general matrix form by setting

where the j ’ s have zero elements except for the following: 1 has diagonal elements 0, r 1 t 1 r , …, r n 1 t 1 r ; 2 has ones on the diagonal above the main diagonal; the last row of 3 consists of the block matrix a n t , a n 1 t ,…, a 1 t ). Here, solutions y and their first n-1 many derivatives are given by the respective components of w as y j 1 t j 1 r w n j 1 : j 1, 2, …, n.

Consider a system of the form w’=t r A w where the elements of A are analytic near infinity and A A 0 t 1 A 1 +…+ t r A r 1 (formal series) where A 0 is a diagonal matrix for distinct complex numbers Λ j along as the diagonal elements. Then, there are constant matrices R, Q j : j 0, 1, …, r 1 and P j : j 0, 1, … so that the columns of (3) t j 0 t j P j t R EXP j 0 r 1 t j Q j are each asymptotic series for some solution of (1). Here the matrices R, Q j are each diagonal, P 0 may be taken as the identity matrix, and Q 0 may be taken to equal A 0 ; here, EXP denotes the matrix exponential and t R EXP ln t R . This is a special case of Theorems 2.1 and 4.1 of Chapt. 5 [1]. (In fact, this method produces exact solutions in cases not treated here. Cf. Chapt. 4 [1].) We will not repeat the entire proof of this result here, but we will elaborat on the construction of the Q j ’s and R. By methods which amount to a matrix version of the dominant balance method, the matrices P j , Q j 1 : 1 j r 1 and R satisfy

For j 1, let us denote by P j are the same matrix as P j ’s above but with zeros on their main diagonals. For our objectives, we need only to solve for the corresponding P j ’s to proceed: We obtain (4)

By ‘Diag( )’ we mean that diagonal matrix whose elements match those of the argument along the main diagonal: This is not exactly the same as the Mathematica algorithm ‘Diagonal[ ]’. Then for each 2 k r (provided r>1) we may recursively compute the Q k ’ s, P k ’ s, and R by ( 5)

We begin with an example involving rational functions of t. (Such examples suffice in our general setting by our hypothesis on A t ). Consider the case r=1, x 0 0, y 0 4, and (8) A t 0.1 0 0 0.5 t 1 1 1 1 0 t 2 0 1 0 1 . We set

In [1]:= A0 1, 0 , 0, 0.5 ; Q0 A0; P0 IdentityMatrix 2 ; A1 1, 1 , 1, 0 ;

as we need only terms of order t and 1 for accuracy up to O t 1 in the first factor of (3).

We follow equations (4) to solve Q 1 and P 1 (we will denote respective P j ’s by Pj’ s in our inputs cells) and ( 7) to solve R: Asymptotic solutions up to first derivatives are given by the columns of the matrix V which we compute by

In [7] We find that there are two linearly independent solution vectors V 1 and V 2 satisfying

We compare our result with exact solution of an initial value problem of the following asymptotically equivalent [2] system:

Routine calculations, involving hypergeometric functions, show that x t = C 0.25 t 2 2 t 1 o t t 2 1 o t for a constant C 12.09.

A natural way to check our results is to compute numerical results via NDSolve[ ].

Here, we compare our dominant asymptotic calculations (from V 2 ) to numerical results (via Interpolation) of the original equation. Below we expect our graphs to have horizontal asymptotes; instead, we find the calculations to be unstable for t 9.

In [8]

We can apply these asymptotic estimates to produce stable solutions from ND-

Solve[ ] for large t. We will replace x t by 2 te .25 t 2 x t and y t by t 2 e .25 t 2 y t to produce a solution x, y with less drastic decay rates: The solutions seem to be tending to limiting values for t near 600 : All with no call to any special calculational options.

It is beyond the scope of this article to compare and contrast numerical schemes since our objective is to predict theoretically the growth behavior of solution as demonstrated using default options of NDSolve[ ]. Our method, we note, may serve to improve stability and reduce stiffness, yet we do not quantify

…(Full text truncated)…