6.13 airy’s equation (7.12.00)

6.13.1 Chuck Baker

Is it possible to develop a power series solution, in powers of x, for Airy’s equations using Maple?

\[ y'' - x y = 0 \]

6.13.2 Chris Eilbeck (8.12.00)

Use the expansions given in Abramowitz and Stegun, Handbook of Mathematical Functions, Section 10.4.

6.13.3 Bill Bauldry (8.12.00)

A naive approach will work here to get one solution:

>  de := (D@@2)(y)(x) - x*y(x) =0; 
 
                            (2) 
                    de := (D   )(y)(x) - x y(x) = 0 
 
>  Order := 10: # if you want more than 5 terms 
>  dsolve(de, y(x), series); 
 
                                       3                 4 
   y(x) = y(0) + D(y)(0) x + 1/6 y(0) x  + 1/12 D(y)(0) x  + 
 
                     6                  7                 9      10 
         1/180 y(0) x  + 1/504 D(y)(0) x  + 1/12960 y(0) x  + O(x  )
 

But it may be more fun to work with the predefined functions:

>  Order := 6: 
>  taylor(AiryAi(x), x); 
>  taylor(AiryBi(x), x); 
 
          (1/3)          (1/6)                        (1/3) 
         3              3      GAMMA(2/3)            3         3 
   1/3 ---------- - 1/2 ----------------- x + 1/18 ---------- x  - 
       GAMMA(2/3)              Pi                  GAMMA(2/3) 
 
               (1/6) 
              3      GAMMA(2/3)  4      6 
         1/24 ----------------- x  + O(x ) 
                     Pi 
 
          (5/6)          (2/3)                        (5/6) 
         3              3      GAMMA(2/3)            3         3 
   1/3 ---------- + 1/2 ----------------- x + 1/18 ---------- x  + 
       GAMMA(2/3)              Pi                  GAMMA(2/3) 
 
               (2/3) 
              3      GAMMA(2/3)  4      6 
         1/24 ----------------- x  + O(x ) 
                     Pi
 

6.13.4 Dr Francis J. Wright (8.12.00)

From the online help for AiryAi in maple6:

series(AiryAi(x),x,4); 
 
         1/3           1/6                         1/3 
        3             3    GAMMA(2/3)             3        3      4 
 1/3 ---------- - 1/2 --------------- x + 1/18 ---------- x  + O(x ) 
     GAMMA(2/3)             Pi                 GAMMA(2/3)
 

There are probably also harder ways to do it!

6.13.5 Craig B. Watkins (8.12.00)

It’s quite easy. One example I have is:

Order:=20; 
f1:=dsolve({diff(y(x),x)=exp(x*y),y(0)=0},y(x),'type=series'); 
f2:=subs(f1,y(x)); 
f3:=convert(f2,polynom); 
plot(f3,x=-1..1); 
g1:=dsolve({diff(y(x),x$2)+x*y=0,y(0)=0,D(y)(0)=1},y(x),'type=series'); 
g2:=subs(g1,y(x)); 
g3:=convert(g2,polynom); 
plot(g3,x=-1..5,scaling=constrained);
 

This may be more than you wanted (and I use \(y'' + x y = 0\) for this example), but it’s what I like to show undergraduates.

An explanation of what this is used for (somewhat MIT-dependent in terms of subject matter is at

http://web.mit.edu/18.03-esg/www/cws00/maple/seriessol.html

6.13.6 Willard, Daniel Dr (8.12.00)

Try this (Release 5.1):

 restart;#?dsolve[series] 
 with(ODEtools); 
 ode:=diff(y(x),x$2)-x*y(x)=0; 
 Order:=20; 
 dsolve(ode,y(x),'type=series'); 
  # or dsolve({ode, y(0)=A,D(y)(0)=B},y(x),'type=series');
 

6.13.7 Robert Israel (8.12.00)

There were lots of replies (including mine) that used the "series" option for dsolve, but Chuck later confirmed to me that he was interested in a symbolic expression for the coefficients. In Maple 6, this can be obtained using the new "formal_series" option:

 
 dsolve((D@@2)(y)(x)-x*y(x)=0,y(x),'formal_series',coeffs=mhypergeom); 
 
                         4 
  {_C[1] x + 1/12 _C[1] x 
 
                                /infinity                               \ 
                                | -----              n  (3 n + 1)       | 
                                |  \            (1/3)  x                | 
           2/9 _C[1] Pi sqrt(3) |   )     ------------------------------| 
                                |  /       n                            | 
                                | -----   3  GAMMA(n + 4/3) GAMMA(n + 1)| 
                                \ n = 2                                 / 
         + --------------------------------------------------------------, 
                                     GAMMA(2/3) 
 
                           3 
        _C[1] + 1/6 _C[1] x 
 
                            /infinity                               \ 
                            | -----                n  (3 n)         | 
                            |  \              (1/3)  x              | 
         + _C[1] GAMMA(2/3) |   )     ------------------------------|} 
                            |  /                    n               | 
                            | -----   GAMMA(n + 1) 3  GAMMA(n + 2/3)| 
                            \ n = 2                                 /
 

6.13.8 Francois Boucher (8.12.00)

Two ideas (with mapleV4) :

> with(powseries): 
> edo:=diff(y(x),x,x)-x*y(x); 
 
                             / 2      \ 
                             |d       | 
                      edo := |--- y(x)| - x y(x) 
                             |  2     | 
                             \dx      / 
 
> sol:=powsolve({edo,y(0)=1,D(y)(0)=0}); 
 
                    sol := proc(powparm)  ...  end 
 
 
 sol(_k); 
 
                              a(_k - 3) 
                             ----------- 
                             _k (_k - 1) 
 
> ak:=rsolve({a(k)=a(k-1)/3/k/(3*k-1),a(0)=1},a(k)); 
 
                              (-k) 
                             9     GAMMA(2/3) 
                  ak := --------------------------- 
                        GAMMA(k + 1) GAMMA(k + 2/3) 
 
> sum(ak*x^(3 *k),k=0..infinity); 
 
                                               3 
                     hypergeom([], [2/3], 1/9 x ) 
 
> restart: 
> with(share): 
#See ?share and ?share,contents for information about the share library 
> readshare(gfun,analysis): 
> with(gfun): 
> edo:=diff(y(x),x,x)-x*y(x); 
 
                             / 2      \ 
                             |d       | 
                      edo := |--- y(x)| - x y(x) 
                             |  2     | 
                             \dx      / 
 
> rec:=diffeqtorec({edo,y(0)=1,D(y)(0)=0},y(x),u(m)); 
 
rec := 
               2 
    {-u(m) + (m  + 5 m + 6) u(m + 3), u(2) = 0, u(0) = 1, u(1) = 0} 
rsolve don't work on this but : 
 
> rec2:=op(1,rec); 
 
                                 2 
               rec2 := -u(m) + (m  + 5 m + 6) u(m + 3) 
 
> rec2:=subs([u(m)=v(n),u(m+3)=v(n+1),m=3*n],rec2); 
 
                                  2 
              rec2 := -v(n) + (9 n  + 15 n + 6) v(n + 1) 
 
> vn:=rsolve({rec2=0,v(0)=1},v(n)); 
 
                              (-n) 
                             9     GAMMA(2/3) 
                  vn := --------------------------- 
                        GAMMA(n + 1) GAMMA(n + 2/3) 
 
> sol:=sum(vn*x^(3*n),n=0..infinity); 
 
                                                  3 
                 sol := hypergeom([], [2/3], 1/9 x )