#set document(title: "7.3 Singular Points and the Method of Frobenius", author: "Jiří Lebl") #set page(width: 8.5in, height: auto, margin: 1in) #import "@preview/cetz:0.5.2" #set text(font: ("STIX Two Text", "Libertinus Serif", "New Computer Modern"), size: 10.5pt, lang: "en") #show math.equation: set text(font: ("STIX Two Math", "New Computer Modern Math")) #set par(justify: true, leading: 0.62em, spacing: 0.9em) #set enum(spacing: 1.1em) // room between list items so tall inline fractions don't collide #set list(spacing: 1.1em) #set table(stroke: 0.5pt + rgb("#c7ccd3")) #let BLUE = rgb("#183B6F") // brand navy — section bars + example/solution labels (white on navy 11.09:1) #let ORANGE = rgb("#A94509") // brand primary-700 — AA-safe deep orange for TEXT (5.93:1 on white; raw brand #F37021 is 2.94:1 and must never carry text) #let RED = rgb("#DC2626") // brand error-600 #let GREEN = rgb("#059669") // brand success-600 (decoration only; small green text uses green-text #007942) #show heading.where(level: 1): it => block(width: 100%, above: 0pt, below: 16pt, fill: gradient.linear(BLUE, rgb("#2C5AA0")), inset: (x: 14pt, y: 12pt), radius: 3pt, text(fill: white, weight: "bold", size: 19pt, it.body)) #show heading.where(level: 2): it => block(width: 100%, above: 18pt, below: 10pt, fill: BLUE, inset: (x: 10pt, y: 6pt), radius: 2pt, text(fill: white, weight: "bold", size: 12pt, it.body)) #show heading.where(level: 3): it => text(fill: ORANGE, weight: "bold", size: 12.5pt, it.body) #show heading.where(level: 4): it => text(fill: BLUE, weight: "bold", size: 10.5pt, it.body) #let examplebox(label, title, body) = block(width: 100%, breakable: true, fill: rgb("#EFF1F5"), stroke: 0.5pt + rgb("#CFDDF0"), radius: 4pt, inset: 10pt, above: 12pt, below: 12pt)[ #block(below: 6pt)[#box(fill: BLUE, inset: (x: 6pt, y: 2pt), radius: 2pt, text(fill: white, weight: "bold", size: 8.5pt, label)) #h(0.4em) #strong[#title]] #body] // rail = decorative left rule (raw brand token); labelcolor = AA-safe label text shade #let notebox(label, rail, labelcolor, tint, body) = block(width: 100%, breakable: true, fill: tint, stroke: (left: 3pt + rail), inset: (left: 10pt, rest: 8pt), radius: (right: 4pt), above: 11pt, below: 11pt)[ #text(fill: labelcolor, weight: "bold", size: 7.5pt, tracking: 0.5pt)[#upper(label)] #linebreak() #body] #let solutionbox(body) = block(above: 4pt, below: 8pt)[ #text(fill: BLUE, weight: "bold", size: 8.5pt)[Solution] #linebreak() #body] #let figph(msg) = block(width: 100%, height: 60pt, fill: rgb("#f6f7f9"), stroke: (paint: rgb("#c7ccd3"), dash: "dashed"), radius: 4pt, inset: 10pt)[ #align(center + horizon, text(fill: rgb("#889"), style: "italic", size: 9pt, msg))] // Standardize inlined figure sizes: measure the natural CeTZ canvas, then scale to a // consistent envelope (aspect-aware; see build_typst.py FIG_* constants). Unlike the // print preamble, dimensions are FLOORED: in an editor a user can trim a figure to a // degenerate 1-D shape (a bare line), and w/h or tw/w would then divide by zero. #let _STD_W = 3.5 #let _WIDE_W = 5.6 #let _MAX_H = 3.4 #let _ASPECT_WIDE = 2.2 #let _UPSCALE_MAX = 1.15 #let stdfig(body) = context { let m = measure(body) let w = calc.max(m.width / 1in, 0.01) let h = calc.max(m.height / 1in, 0.01) let tw = if w / h > _ASPECT_WIDE { _WIDE_W } else { _STD_W } let s = calc.min(tw / w, _MAX_H / h, _UPSCALE_MAX) align(center, box(scale(x: s * 100%, y: s * 100%, reflow: true, body))) } #show figure: set block(breakable: false) #set figure(gap: 8pt) #show figure.caption: set text(size: 8.5pt, fill: rgb("#555")) == 7.3#h(0.6em)Singular Points and the Method of Frobenius === Examples While behavior of ODEs at singular points is more complicated, certain singular points are not especially difficult to solve. Let us look at some examples before giving a general method. We may be lucky and obtain a power series solution using the method of the previous section, but in general we may have to try other things. The slope field of x y′ = r y with three of its solutions y = C|x|^r drawn through it. The equation worked just below, 2xy′ − y = 0, is the case r = 1/2 — where the figure opens — because substituting y = x^r turns the ODE into the indicial equation for r. Drag the slider to watch the singular point change character: for r \< 0 the solutions themselves blow up at x = 0; at r = 0 they flatten into the horizontal lines y = C; for 0 \< r \< 1, the range containing Lebl's r = 1/2, they reach the origin with a vertical tangent, which is why the trivial y = 0 is the only solution differentiable there; at r = 1 the two branches meet in a corner; and only for r \> 1 do they arrive with zero slope and the singular point stop hurting. {"functions":\[{"color":"\#1f4e79","domain":null,"expression":"r\*y/x","expression2":"","id":"frob-field","inequality":"lt","lineDash":\[\],"lineWidth":2,"mode":"slope-field","tMax":10,"tMin":-10,"variable":"x","visible":true},{"color":"\#b8442a","domain":null,"expression":"abs(x)^r","expression2":"","id":"frob-sol-1","inequality":"lt","lineDash":\[\],"lineWidth":2.5,"mode":"cartesian","tMax":10,"tMin":-10,"variable":"x","visible":true},{"color":"\#2e6b3f","domain":null,"expression":"0.5\*abs(x)^r","expression2":"","id":"frob-sol-half","inequality":"lt","lineDash":\[\],"lineWidth":2.5,"mode":"cartesian","tMax":10,"tMin":-10,"variable":"x","visible":true},{"color":"\#6b3fa0","domain":null,"expression":"-0.8\*abs(x)^r","expression2":"","id":"frob-sol-neg","inequality":"lt","lineDash":\[\],"lineWidth":2.5,"mode":"cartesian","tMax":10,"tMin":-10,"variable":"x","visible":true}\],"grid":{"axisColor":"\#333333","gridColor":"\#cccccc","majorSpacing":1,"minorGridColor":"\#eeeeee","showAxes":true,"showLabels":true},"parameters":\[{"id":"frob-r","label":"Frobenius exponent r","max":2,"min":-2,"name":"r","step":0.1,"unit":"","value":0.5}\],"title":"Slope field for x y' = r y at a singular point","version":1,"viewport":{"centerX":0,"centerY":0,"scale":95}} #examplebox("Example 1")[][ Let us first look at a simple first order equation #math.equation(block: true, alt: "2 x y prime minus y equals 0 .")[$2 x y^(′) − y = 0 .$] Note that #math.equation(block: false, alt: "x equals 0")[$x = 0$] is a singular point. If we only try to plug in #math.equation(block: true, alt: "y equals ∑ k equals 0 infinity a sub k x to the power k ,")[$y = ∑_(k = 0)^(∞) a_(k) x^(k) ,$] we obtain #math.equation(block: true, alt: "0 equals 2 x y prime minus y equals 2 x open parenthesis ∑ k equals 1 infinity k a sub k x to the power k minus 1 close parenthesis minus open parenthesis ∑ k equals 0 infinity a sub k x to the power k close parenthesis; equals a sub 0 plus ∑ k equals 1 infinity open parenthesis 2 k a sub k minus a sub k close parenthesis x to the power k .")[$0 = 2 x y^(′) − y = 2 x ( ∑_(k = 1)^(∞) k a_(k) x^(k − 1) ) − ( ∑_(k = 0)^(∞) a_(k) x^(k) ) \ = a_(0) + ∑_(k = 1)^(∞) ( 2 k a_(k) − a_(k) ) x^(k) .$] First, #math.equation(block: false, alt: "a sub 0 equals 0")[$a_(0) = 0$]. Next, the only way to solve #math.equation(block: false, alt: "0 equals 2 k a sub k minus a sub k equals open parenthesis 2 k minus 1 close parenthesis a sub k")[$0 = 2 k a_(k) − a_(k) = ( 2 k − 1 ) a_(k)$] for #math.equation(block: false, alt: "k equals 1 , 2 , 3 , and so on")[$k = 1 , 2 , 3 , …$] is for #math.equation(block: false, alt: "a sub k equals 0")[$a_(k) = 0$] for all #math.equation(block: false, alt: "k")[$k$]. Therefore we only get the trivial solution #math.equation(block: false, alt: "y equals 0")[$y = 0$]. We need a nonzero solution to get the general solution. Let us try #math.equation(block: false, alt: "y equals x to the power r")[$y = x^(r)$] for some real number #math.equation(block: false, alt: "r")[$r$]. Consequently our solution---if we can find one---may only make sense for positive #math.equation(block: false, alt: "x")[$x$]. Then #math.equation(block: false, alt: "y prime equals r x to the power r minus 1")[$y^(′) = r x^(r − 1)$]. So #math.equation(block: true, alt: "0 equals 2 x y prime minus y equals 2 x r x to the power r minus 1 minus x to the power r equals open parenthesis 2 r minus 1 close parenthesis x to the power r .")[$0 = 2 x y^(′) − y = 2 x r x^(r − 1) − x^(r) = ( 2 r − 1 ) x^(r) .$] Therefore #math.equation(block: false, alt: "r equals the fraction 1 over 2")[$r = frac(1, 2)$], or in other words #math.equation(block: false, alt: "y equals x to the power 1 / 2")[$y = x^(1 / 2)$]. Multiplying by a constant, the general solution for positive #math.equation(block: false, alt: "x")[$x$] is #math.equation(block: true, alt: "y equals C x to the power 1 / 2 .")[$y = C x^(1 / 2) .$] If #math.equation(block: false, alt: "C ⧸ equals 0")[$C "⧸" = 0$] then the derivative of the solution "blows up" at #math.equation(block: false, alt: "x equals 0")[$x = 0$] (the singular point). There is only one solution that is differentiable at #math.equation(block: false, alt: "x equals 0")[$x = 0$] and that's the trivial solution #math.equation(block: false, alt: "y equals 0")[$y = 0$]. ] Not every problem with a singular point has a solution of the form #math.equation(block: false, alt: "y equals x to the power r")[$y = x^(r)$], of course. But perhaps we can combine the methods. What we will do is to try a solution of the form #math.equation(block: true, alt: "y equals x to the power r f open parenthesis x close parenthesis")[$y = x^(r) f ( x )$] where #math.equation(block: false, alt: "f open parenthesis x close parenthesis")[$f ( x )$] is an analytic function. #examplebox("Example 2")[][ Suppose that we have the equation #math.equation(block: true, alt: "4 x squared y double prime minus 4 x squared y prime plus open parenthesis 1 minus 2 x close parenthesis y equals 0 ,")[$4 x^(2) y^(″) − 4 x^(2) y^(′) + ( 1 − 2 x ) y = 0 ,$] and again note that #math.equation(block: false, alt: "x equals 0")[$x = 0$] is a singular point. Let us try #math.equation(block: true, alt: "y equals x to the power r ∑ k equals 0 infinity a sub k x to the power k equals ∑ k equals 0 infinity a sub k x to the power k plus r ,")[$y = x^(r) ∑_(k = 0)^(∞) a_(k) x^(k) = ∑_(k = 0)^(∞) a_(k) x^(k + r) ,$] where #math.equation(block: false, alt: "r")[$r$] is a real number, not necessarily an integer. Again if such a solution exists, it may only exist for positive #math.equation(block: false, alt: "x")[$x$]. First let us find the derivatives #math.equation(block: true, alt: "y prime equals ∑ k equals 0 infinity open parenthesis k plus r close parenthesis a sub k x to the power k plus r minus 1 ,; y double prime equals ∑ k equals 0 infinity open parenthesis k plus r close parenthesis open parenthesis k plus r minus 1 close parenthesis a sub k x to the power k plus r minus 2 .")[$y^(′) = ∑_(k = 0)^(∞) ( k + r ) a_(k) x^(k + r − 1) , \ y^(″) = ∑_(k = 0)^(∞) ( k + r ) ( k + r − 1 ) a_(k) x^(k + r − 2) .$] Plugging Equations \\ref{ex2eq2} - \\ref{ex2eq3} into our original differential equation (Equation \\ref{ex2eq1}) we obtain #math.equation(block: true, alt: "0 equals 4 x squared y double prime minus 4 x squared y prime plus open parenthesis 1 minus 2 x close parenthesis y; equals 4 x squared open parenthesis ∑ k equals 0 infinity open parenthesis k plus r close parenthesis open parenthesis k plus r minus 1 close parenthesis a sub k x to the power k plus r minus 2 close parenthesis minus 4 x squared open parenthesis ∑ k equals 0 infinity open parenthesis k plus r close parenthesis a sub k x to the power k plus r minus 1 close parenthesis plus open parenthesis 1 minus 2 x close parenthesis open parenthesis ∑ k equals 0 infinity a sub k x to the power k plus r close parenthesis; equals open parenthesis ∑ k equals 0 infinity 4 open parenthesis k plus r close parenthesis open parenthesis k plus r minus 1 close parenthesis a sub k x to the power k p")[$0 = 4 x^(2) y^(″) − 4 x^(2) y^(′) + ( 1 − 2 x ) y \ = 4 x^(2) ( ∑_(k = 0)^(∞) ( k + r ) ( k + r − 1 ) a_(k) x^(k + r − 2) ) − 4 x^(2) ( ∑_(k = 0)^(∞) ( k + r ) a_(k) x^(k + r − 1) ) + ( 1 − 2 x ) ( ∑_(k = 0)^(∞) a_(k) x^(k + r) ) \ = ( ∑_(k = 0)^(∞) 4 ( k + r ) ( k + r − 1 ) a_(k) x^(k + r) ) − ( ∑_(k = 0)^(∞) 4 ( k + r ) a_(k) x^(k + r + 1) ) + ( ∑_(k = 0)^(∞) a_(k) x^(k + r) ) − ( ∑_(k = 0)^(∞) 2 a_(k) x^(k + r + 1) ) \ = ( ∑_(k = 0)^(∞) 4 ( k + r ) ( k + r − 1 ) a_(k) x^(k + r) ) − ( ∑_(k = 1)^(∞) 4 ( k + r − 1 ) a_(k − 1) x^(k + r) ) + ( ∑_(k = 0)^(∞) a_(k) x^(k + r) ) − ( ∑_(k = 1)^(∞) 2 a_(k − 1) x^(k + r) ) \ = 4 r ( r − 1 ) a_(0) x^(r) + a_(0) x^(r) + ∑_(k = 1)^(∞) ( 4 ( k + r ) ( k + r − 1 ) a_(k) − 4 ( k + r − 1 ) a_(k − 1) + a_(k) − 2 a_(k − 1) ) x^(k + r) \ = ( 4 r ( r − 1 ) + 1 ) a_(0) x^(r) + ∑_(k = 1)^(∞) ( ( 4 ( k + r ) ( k + r − 1 ) + 1 ) a_(k) − ( 4 ( k + r − 1 ) + 2 ) a_(k − 1) ) x^(k + r) .$] To have a solution we must first have #math.equation(block: false, alt: "open parenthesis 4 r open parenthesis r minus 1 close parenthesis plus 1 close parenthesis a sub 0 equals 0")[$( 4 r ( r − 1 ) + 1 ) a_(0) = 0$]. Supposing that #math.equation(block: false, alt: "a sub 0 ⧸ equals 0")[$a_(0) "⧸" = 0$] we obtain #math.equation(block: true, alt: "4 r open parenthesis r minus 1 close parenthesis plus 1 equals 0 .")[$4 r ( r − 1 ) + 1 = 0 .$] This equation is called the #emph[indicial equation]. This particular indicial equation has a double root at #math.equation(block: false, alt: "r equals the fraction 1 over 2")[$r = frac(1, 2)$]. OK, so we know what #math.equation(block: false, alt: "r")[$r$] has to be. That knowledge we obtained simply by looking at the coefficient of #math.equation(block: false, alt: "x to the power r")[$x^(r)$]. All other coefficients of #math.equation(block: false, alt: "x to the power k plus r")[$x^(k + r)$] also have to be zero so #math.equation(block: true, alt: "open parenthesis 4 open parenthesis k plus r close parenthesis open parenthesis k plus r minus 1 close parenthesis plus 1 close parenthesis a sub k minus open parenthesis 4 open parenthesis k plus r minus 1 close parenthesis plus 2 close parenthesis a sub k minus 1 equals 0 .")[$( 4 ( k + r ) ( k + r − 1 ) + 1 ) a_(k) − ( 4 ( k + r − 1 ) + 2 ) a_(k − 1) = 0 .$] If we plug in #math.equation(block: false, alt: "r equals the fraction 1 over 2")[$r = frac(1, 2)$] and solve for #math.equation(block: false, alt: "a sub k")[$a_(k)$] we get #math.equation(block: true, alt: "a sub k equals the fraction 4 open parenthesis k plus the fraction 1 over 2 minus 1 close parenthesis plus 2 over 4 open parenthesis k plus the fraction 1 over 2 close parenthesis open parenthesis k plus the fraction 1 over 2 minus 1 close parenthesis plus 1 a sub k minus 1 equals the fraction 1 over k a sub k minus 1 .")[$a_(k) = frac(4 ( k + frac(1, 2) − 1 ) + 2, 4 ( k + frac(1, 2) ) ( k + frac(1, 2) − 1 ) + 1) a_(k − 1) = frac(1, k) a_(k − 1) .$] Let us set #math.equation(block: false, alt: "a sub 0 equals 1")[$a_(0) = 1$]. Then #math.equation(block: true, alt: "a sub 1 equals the fraction 1 over 1 a sub 0 equals 1 , a sub 2 equals the fraction 1 over 2 a sub 1 equals the fraction 1 over 2 , a sub 3 equals the fraction 1 over 3 a sub 2 equals the fraction 1 over 3 times 2 , a sub 4 equals the fraction 1 over 4 a sub 3 equals the fraction 1 over 4 times 3 times 2 , and so on")[$a_(1) = frac(1, 1) a_(0) = 1 , #h(2em) a_(2) = frac(1, 2) a_(1) = frac(1, 2) , #h(2em) a_(3) = frac(1, 3) a_(2) = frac(1, 3 · 2) , #h(2em) a_(4) = frac(1, 4) a_(3) = frac(1, 4 · 3 · 2) , #h(2em) …$] Extrapolating, we notice that #math.equation(block: true, alt: "a sub k equals the fraction 1 over k open parenthesis k minus 1 close parenthesis open parenthesis k minus 2 close parenthesis ⋯ 3 times 2 equals the fraction 1 over k ! .")[$a_(k) = frac(1, k ( k − 1 ) ( k − 2 ) ⋯ 3 · 2) = frac(1, k !) .$] In other words, #math.equation(block: true, alt: "y equals ∑ k equals 0 infinity a sub k x to the power k plus r equals ∑ k equals 0 infinity the fraction 1 over k ! x to the power k plus 1 / 2 equals x to the power 1 / 2 ∑ k equals 0 infinity the fraction 1 over k ! x to the power k equals x to the power 1 / 2 e to the power x .")[$y = ∑_(k = 0)^(∞) a_(k) x^(k + r) = ∑_(k = 0)^(∞) frac(1, k !) x^(k + 1 / 2) = x^(1 / 2) ∑_(k = 0)^(∞) frac(1, k !) x^(k) = x^(1 / 2) e^(x) .$] That was lucky! In general, we will not be able to write the series in terms of elementary functions. We have one solution, let us call it #math.equation(block: false, alt: "y sub 1 equals x to the power 1 / 2 e to the power x")[$y_(1) = x^(1 / 2) e^(x)$]. But what about a second solution? If we want a general solution, we need two linearly independent solutions. Picking #math.equation(block: false, alt: "a sub 0")[$a_(0)$] to be a different constant only gets us a constant multiple of #math.equation(block: false, alt: "y sub 1")[$y_(1)$], and we do not have any other #math.equation(block: false, alt: "r")[$r$] to try; we only have one solution to the indicial equation. Well, there are powers of #math.equation(block: false, alt: "x")[$x$] floating around and we are taking derivatives, perhaps the logarithm (the antiderivative of #math.equation(block: false, alt: "x to the power minus 1")[$x^(− 1)$]) is around as well. It turns out we want to try for another solution of the form #math.equation(block: true, alt: "y sub 2 equals ∑ k equals 0 infinity b sub k x to the power k plus r plus open parenthesis ln x close parenthesis y sub 1 ,")[$y_(2) = ∑_(k = 0)^(∞) b_(k) x^(k + r) + ( ln x ) y_(1) ,$] which in our case is #math.equation(block: true, alt: "y sub 2 equals ∑ k equals 0 infinity b sub k x to the power k plus 1 / 2 plus open parenthesis ln x close parenthesis x to the power 1 / 2 e to the power x .")[$y_(2) = ∑_(k = 0)^(∞) b_(k) x^(k + 1 / 2) + ( ln x ) x^(1 / 2) e^(x) .$] We now differentiate this equation, substitute into the differential equation and solve for #math.equation(block: false, alt: "b sub k")[$b_(k)$]. A long computation ensues and we obtain some recursion relation for #math.equation(block: false, alt: "b sub k")[$b_(k)$]. The reader can (and should) try this to obtain for example the first three terms #math.equation(block: true, alt: "b sub 1 equals b sub 0 minus 1 , b sub 2 equals the fraction 2 b sub 1 minus 1 over 4 , b sub 3 equals the fraction 6 b sub 2 minus 1 over 18 , …")[$b_(1) = b_(0) − 1 , #h(2em) b_(2) = frac(2 b_(1) − 1, 4) , #h(2em) b_(3) = frac(6 b_(2) − 1, 18) , #h(2em) …$] We then fix #math.equation(block: false, alt: "b sub 0")[$b_(0)$] and obtain a solution #math.equation(block: false, alt: "y sub 2")[$y_(2)$]. Then we write the general solution as #math.equation(block: false, alt: "y equals A y sub 1 plus B y sub 2")[$y = A y_(1) + B y_(2)$]. ] === Method of Frobenius Before giving the general method, let us clarify when the method applies. Let #math.equation(block: true, alt: "p open parenthesis x close parenthesis y double prime plus q open parenthesis x close parenthesis y prime plus r open parenthesis x close parenthesis y equals 0")[$p ( x ) y^(″) + q ( x ) y^(′) + r ( x ) y = 0$] be an ODE. As before, if #math.equation(block: false, alt: "p open parenthesis x sub 0 close parenthesis equals 0")[$p ( x_(0) ) = 0$], then #math.equation(block: false, alt: "x sub 0")[$x_(0)$] is a singular point. If, furthermore, the limits #math.equation(block: true, alt: "lim sub x → x sub 0 open parenthesis x minus x sub 0 close parenthesis the fraction q open parenthesis x close parenthesis over p open parenthesis x close parenthesis and lim sub x → x sub 0 open parenthesis x minus x sub 0 close parenthesis squared the fraction r open parenthesis x close parenthesis over p open parenthesis x close parenthesis")[$lim_(x → x_(0)) " " ( x − x_(0) ) frac(q ( x ), p ( x )) #h(2em) "and" #h(2em) lim_(x → x_(0)) " " ( x − x_(0) )^(2) frac(r ( x ), p ( x ))$] both exist and are finite, then we say that #math.equation(block: false, alt: "x sub 0")[$x_(0)$] is a #emph[regular singular point]. #examplebox("Example 3")[Expansion around a regular singular point][ Often, and for the rest of this section, #math.equation(block: false, alt: "x sub 0 equals 0")[$x_(0) = 0$]. Consider #math.equation(block: true, alt: "x squared y double prime plus x open parenthesis 1 plus x close parenthesis y prime plus open parenthesis π plus x squared close parenthesis y equals 0 .")[$x^(2) y^(″) + x ( 1 + x ) y^(′) + ( π + x^(2) ) y = 0 .$] Write #math.equation(block: true, alt: "lim sub x → 0 x the fraction q open parenthesis x close parenthesis over p open parenthesis x close parenthesis equals lim sub x → 0 x the fraction x open parenthesis 1 plus x close parenthesis over x squared equals lim sub x → 0 open parenthesis 1 plus x close parenthesis equals 1 ,; lim sub x → 0 x squared the fraction r open parenthesis x close parenthesis over p open parenthesis x close parenthesis equals lim sub x → 0 x squared the fraction open parenthesis π plus x squared close parenthesis over x squared equals lim sub x → 0 open parenthesis π plus x squared close parenthesis equals π")[$lim_(x → 0) " " x frac(q ( x ), p ( x )) = lim_(x → 0) " " x frac(x ( 1 + x ), x^(2)) = lim_(x → 0) " " ( 1 + x ) = 1 , \ lim_(x → 0) " " x^(2) frac(r ( x ), p ( x )) = lim_(x → 0) " " x^(2) frac(( π + x^(2) ), x^(2)) = lim_(x → 0) " " ( π + x^(2) ) = π$] So #math.equation(block: false, alt: "x equals 0")[$x = 0$] is a regular singular point. On the other hand if we make the slight change #math.equation(block: true, alt: "x squared y double prime plus open parenthesis 1 plus x close parenthesis y prime plus open parenthesis π plus x squared close parenthesis y equals 0 ,")[$x^(2) y^(″) + ( 1 + x ) y^(′) + ( π + x^(2) ) y = 0 ,$] then #math.equation(block: true, alt: "lim sub x → 0 x the fraction q open parenthesis x close parenthesis over p open parenthesis x close parenthesis equals lim sub x → 0 x the fraction open parenthesis 1 plus x close parenthesis over x squared equals lim sub x → 0 the fraction 1 plus x over x equals DNE .")[$lim_(x → 0) " " x frac(q ( x ), p ( x )) = lim_(x → 0) " " x frac(( 1 + x ), x^(2)) = lim_(x → 0) " " frac(1 + x, x) = "DNE" .$] Here DNE stands for #emph[does not exist]. The point #math.equation(block: false, alt: "0")[$0$] is a singular point, but not a regular singular point. ] Let us now discuss the general #emph[Method of Frobenius]#math.equation(block: false, alt: "to the power 1")[$1$]. Let us only consider the method at the point #math.equation(block: false, alt: "x equals 0")[$x = 0$] for simplicity. The main idea is the following theorem. #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[ #strong[Method of Frobenius] Suppose that #math.equation(block: true, alt: "p open parenthesis x close parenthesis y double prime plus q open parenthesis x close parenthesis y prime plus r open parenthesis x close parenthesis y equals 0")[$p ( x ) y^(″) + q ( x ) y^(′) + r ( x ) y = 0$] has a regular singular point at #math.equation(block: false, alt: "x equals 0")[$x = 0$], then there exists at least one solution of the form #math.equation(block: true, alt: "y equals x to the power r ∑ k equals 0 infinity a sub k x to the power k .")[$y = x^(r) ∑_(k = 0)^(∞) a_(k) x^(k) .$] A solution of this form is called a #emph[Frobenius-type solution]. ] The method usually breaks down like this. + We seek a Frobenius-type solution of the form #math.equation(block: true, alt: "y equals ∑ k equals 0 infinity a sub k x to the power k plus r .")[$y = ∑_(k = 0)^(∞) a_(k) x^(k + r) .$] We plug this #math.equation(block: false, alt: "y")[$y$] into equation (7.3.7). We collect terms and write everything as a single series. + The obtained series must be zero. Setting the first coefficient (usually the coefficient of #math.equation(block: false, alt: "x to the power r")[$x^(r)$]) in the series to zero we obtain the #emph[indicial equation], which is a quadratic polynomial in #math.equation(block: false, alt: "r")[$r$]. + If the indicial equation has two real roots #math.equation(block: false, alt: "r sub 1")[$r_(1)$] and #math.equation(block: false, alt: "r sub 2")[$r_(2)$] such that #math.equation(block: false, alt: "r sub 1 minus r sub 2")[$r_(1) − r_(2)$] is not an integer, then we have two linearly independent Frobenius-type solutions. Using the first root, we plug in #math.equation(block: true, alt: "y sub 1 equals x to the power r sub 1 ∑ k equals 0 infinity a sub k x to the power k ,")[$y_(1) = x^(r_(1)) ∑_(k = 0)^(∞) a_(k) x^(k) ,$] and we solve for all #math.equation(block: false, alt: "a sub k")[$a_(k)$] to obtain the first solution. Then using the second root, we plug in #math.equation(block: true, alt: "y sub 2 equals x to the power r sub 2 ∑ k equals 0 infinity b sub k x to the power k ,")[$y_(2) = x^(r_(2)) ∑_(k = 0)^(∞) b_(k) x^(k) ,$] and solve for all #math.equation(block: false, alt: "b sub k")[$b_(k)$] to obtain the second solution. + If the indicial equation has a doubled root #math.equation(block: false, alt: "r")[$r$], then there we find one solution #math.equation(block: true, alt: "y sub 1 equals x to the power r ∑ k equals 0 infinity a sub k x to the power k ,")[$y_(1) = x^(r) ∑_(k = 0)^(∞) a_(k) x^(k) ,$] and then we obtain a new solution by plugging #math.equation(block: true, alt: "y sub 2 equals x to the power r ∑ k equals 0 infinity b sub k x to the power k plus open parenthesis ln x close parenthesis y sub 1 ,")[$y_(2) = x^(r) ∑_(k = 0)^(∞) b_(k) x^(k) + ( ln x ) y_(1) ,$] into Equation (7.3.7) and solving for the constants #math.equation(block: false, alt: "b sub k")[$b_(k)$]. + If the indicial equation has two real roots such that #math.equation(block: false, alt: "r sub 1 minus r sub 2")[$r_(1) − r_(2)$] is an integer, then one solution is #math.equation(block: true, alt: "y sub 1 equals x to the power r sub 1 ∑ k equals 0 infinity a sub k x to the power k ,")[$y_(1) = x^(r_(1)) ∑_(k = 0)^(∞) a_(k) x^(k) ,$] and the second linearly independent solution is of the form #math.equation(block: true, alt: "y sub 2 equals x to the power r sub 2 ∑ k equals 0 infinity b sub k x to the power k plus C open parenthesis ln x close parenthesis y sub 1 ,")[$y_(2) = x^(r_(2)) ∑_(k = 0)^(∞) b_(k) x^(k) + C ( ln x ) y_(1) ,$] where we plug #math.equation(block: false, alt: "y sub 2")[$y_(2)$] into (7.3.7) and solve for the constants #math.equation(block: false, alt: "b sub k")[$b_(k)$] and #math.equation(block: false, alt: "C")[$C$]. + Finally, if the indicial equation has complex roots, then solving for #math.equation(block: false, alt: "a sub k")[$a_(k)$] in the solution #math.equation(block: true, alt: "y equals x to the power r sub 1 ∑ k equals 0 infinity a sub k x to the power k")[$y = x^(r_(1)) ∑_(k = 0)^(∞) a_(k) x^(k)$] results in a complex-valued function---all the #math.equation(block: false, alt: "a sub k")[$a_(k)$] are complex numbers. We obtain our two linearly independent solutions#math.equation(block: false, alt: "squared")[$2$] by taking the real and imaginary parts of #math.equation(block: false, alt: "y")[$y$]. The main idea is to find at least one Frobenius-type solution. If we are lucky and find two, we are done. If we only get one, we either use the ideas above or even a different method such as reduction of order (Exercise 2.1.8) to obtain a second solution. === Bessel Functions An important class of functions that arises commonly in physics are the #emph[Bessel functions]#math.equation(block: false, alt: "cubed")[$3$]. For example, these functions appear when solving the wave equation in two and three dimensions. First we have #emph[Bessel's equation] of order #math.equation(block: false, alt: "p")[$p$]: #math.equation(block: true, alt: "x squared y double prime plus x y prime plus open parenthesis x squared minus p squared close parenthesis y equals 0 .")[$x^(2) y^(″) + x y^(′) + ( x^(2) − p^(2) ) y = 0 .$] We allow #math.equation(block: false, alt: "p")[$p$] to be any number, not just an integer, although integers and multiples of #math.equation(block: false, alt: "the fraction 1 over 2")[$frac(1, 2)$] are most important in applications. When we plug #math.equation(block: true, alt: "y equals ∑ k equals 0 infinity a sub k x to the power k plus r")[$y = ∑_(k = 0)^(∞) a_(k) x^(k + r)$] into Bessel's equation of order #math.equation(block: false, alt: "p")[$p$] we obtain the indicial equation #math.equation(block: true, alt: "r open parenthesis r minus 1 close parenthesis plus r minus p squared equals open parenthesis r minus p close parenthesis open parenthesis r plus p close parenthesis equals 0 .")[$r ( r − 1 ) + r − p^(2) = ( r − p ) ( r + p ) = 0 .$] Therefore we obtain two roots #math.equation(block: false, alt: "r sub 1 equals p")[$r_(1) = p$] and #math.equation(block: false, alt: "r sub 2 equals minus p")[$r_(2) = − p$]. If #math.equation(block: false, alt: "p")[$p$] is not an integer following the method of Frobenius and setting #math.equation(block: false, alt: "a sub 0 equals 1")[$a_(0) = 1$], we obtain linearly independent solutions of the form #math.equation(block: true, alt: "y sub 1 equals x to the power p ∑ k equals 0 infinity the fraction open parenthesis minus 1 close parenthesis to the power k x to the power 2 k over 2 to the power 2 k k ! open parenthesis k plus p close parenthesis open parenthesis k minus 1 plus p close parenthesis ⋯ open parenthesis 2 plus p close parenthesis open parenthesis 1 plus p close parenthesis ,; y sub 2 equals x to the power minus p ∑ k equals 0 infinity the fraction open parenthesis minus 1 close parenthesis to the power k x to the power 2 k over 2 to the power 2 k k ! open parenthesis k minus p close parenthesis open parenthesis k minus 1 minus p close parenthesis ⋯ open parenthesis 2 minus p close parenthesis open parenthesis 1 minus p close parenthesis .")[$y_(1) = x^(p) ∑_(k = 0)^(∞) frac(( − 1 )^(k) x^(2 k), 2^(2 k) k ! ( k + p ) ( k − 1 + p ) ⋯ ( 2 + p ) ( 1 + p )) , \ y_(2) = x^(− p) ∑_(k = 0)^(∞) frac(( − 1 )^(k) x^(2 k), 2^(2 k) k ! ( k − p ) ( k − 1 − p ) ⋯ ( 2 − p ) ( 1 − p )) .$] #notebox("Your Turn", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[ + Verify that the indicial equation of Bessel's equation of order #math.equation(block: false, alt: "p")[$p$] is #math.equation(block: false, alt: "open parenthesis r minus p close parenthesis open parenthesis r plus p close parenthesis equals 0")[$( r − p ) ( r + p ) = 0$]. + Suppose that #math.equation(block: false, alt: "p")[$p$] is not an integer. Carry out the computation to obtain the solutions #math.equation(block: false, alt: "y sub 1")[$y_(1)$] and #math.equation(block: false, alt: "y sub 2")[$y_(2)$] above. ] Bessel functions will be convenient constant multiples of #math.equation(block: false, alt: "y sub 1")[$y_(1)$] and #math.equation(block: false, alt: "y sub 2")[$y_(2)$]. First we must define the #emph[gamma function] #math.equation(block: true, alt: "Γ open parenthesis x close parenthesis equals ∫ 0 infinity t to the power x minus 1 e to the power minus t d t .")[$Γ ( x ) = ∫_(0)^(∞) t^(x − 1) e^(− t) d t .$] Notice that #math.equation(block: false, alt: "Γ open parenthesis 1 close parenthesis equals 1")[$Γ ( 1 ) = 1$]. The gamma function also has a wonderful property #math.equation(block: true, alt: "Γ open parenthesis x plus 1 close parenthesis equals x Γ open parenthesis x close parenthesis .")[$Γ ( x + 1 ) = x Γ ( x ) .$] From this property, one can show that #math.equation(block: false, alt: "Γ open parenthesis n close parenthesis equals open parenthesis n minus 1 close parenthesis !")[$Γ ( n ) = ( n − 1 ) !$] when #math.equation(block: false, alt: "n")[$n$] is an integer, so the gamma function is a continuous version of the factorial. We compute: #math.equation(block: true, alt: "Γ open parenthesis k plus p plus 1 close parenthesis equals open parenthesis k plus p close parenthesis open parenthesis k minus 1 plus p close parenthesis ⋯ open parenthesis 2 plus p close parenthesis open parenthesis 1 plus p close parenthesis Γ open parenthesis 1 plus p close parenthesis ,; Γ open parenthesis k minus p plus 1 close parenthesis equals open parenthesis k minus p close parenthesis open parenthesis k minus 1 minus p close parenthesis ⋯ open parenthesis 2 minus p close parenthesis open parenthesis 1 minus p close parenthesis Γ open parenthesis 1 minus p close parenthesis .")[$Γ ( k + p + 1 ) = ( k + p ) ( k − 1 + p ) ⋯ ( 2 + p ) ( 1 + p ) Γ ( 1 + p ) , \ Γ ( k − p + 1 ) = ( k − p ) ( k − 1 − p ) ⋯ ( 2 − p ) ( 1 − p ) Γ ( 1 − p ) .$] #notebox("Your Turn", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[ Verify the above identities using #math.equation(block: false, alt: "Γ open parenthesis x plus 1 close parenthesis equals x Γ open parenthesis x close parenthesis")[$Γ ( x + 1 ) = x Γ ( x )$]. ] We define the #emph[Bessel functions of the first kind] of order #math.equation(block: false, alt: "p")[$p$] and #math.equation(block: false, alt: "minus p")[$− p$] as #math.equation(block: true, alt: "J sub p open parenthesis x close parenthesis equals the fraction 1 over 2 to the power p Γ open parenthesis 1 plus p close parenthesis y sub 1 equals ∑ k equals 0 infinity the fraction open parenthesis minus 1 close parenthesis to the power k over k ! Γ open parenthesis k plus p plus 1 close parenthesis open parenthesis the fraction x over 2 close parenthesis to the power 2 k plus p ,; J sub minus p open parenthesis x close parenthesis equals the fraction 1 over 2 to the power minus Γ open parenthesis 1 minus p close parenthesis y sub 2 equals ∑ k equals 0 infinity the fraction open parenthesis minus 1 close parenthesis to the power k over k ! Γ open parenthesis k minus p plus 1 close parenthesis open parenthesis the fraction x over 2 close parenthesis to the power 2 k minus p .")[$J_(p) ( x ) = frac(1, 2^(p) Γ ( 1 + p )) y_(1) = ∑_(k = 0)^(∞) frac(attach(( − 1 ), t: k), k ! Γ ( k + p + 1 )) attach(( frac(x, 2) ), t: 2 k + p) , \ J_(− p) ( x ) = frac(1, 2^(−) Γ ( 1 − p )) y_(2) = ∑_(k = 0)^(∞) frac(attach(( − 1 ), t: k), k ! Γ ( k − p + 1 )) attach(( frac(x, 2) ), t: 2 k − p) .$] As these are constant multiples of the solutions we found above, these are both solutions to Bessel's equation of order #math.equation(block: false, alt: "p")[$p$]. The constants are picked for convenience. When #math.equation(block: false, alt: "p")[$p$] is not an integer, #math.equation(block: false, alt: "J sub p")[$J_(p)$] and #math.equation(block: false, alt: "J sub minus p")[$J_(− p)$] are linearly independent. When #math.equation(block: false, alt: "n")[$n$] is an integer we obtain #math.equation(block: true, alt: "J sub n open parenthesis x close parenthesis equals ∑ k equals 0 infinity the fraction open parenthesis minus 1 close parenthesis to the power k over k ! open parenthesis k plus n close parenthesis ! open parenthesis the fraction x over 2 close parenthesis to the power 2 k plus n .")[$J_(n) ( x ) = ∑_(k = 0)^(∞) frac(attach(( − 1 ), t: k), k ! ( k + n ) !) attach(( frac(x, 2) ), t: 2 k + n) .$] In this case it turns out that #math.equation(block: true, alt: "J sub n open parenthesis x close parenthesis equals open parenthesis minus 1 close parenthesis to the power n J sub minus n open parenthesis x close parenthesis ,")[$J_(n) ( x ) = attach(( − 1 ), t: n) J_(− n) ( x ) ,$] and so we do not obtain a second linearly independent solution. The other solution is the so-called #emph[Bessel function of second kind]. These make sense only for integer orders #math.equation(block: false, alt: "n")[$n$] and are defined as limits of linear combinations of #math.equation(block: false, alt: "J sub p open parenthesis x close parenthesis")[$J_(p) ( x )$] and #math.equation(block: false, alt: "J sub minus p open parenthesis x close parenthesis")[$J_(− p) ( x )$] as #math.equation(block: false, alt: "p")[$p$] approaches #math.equation(block: false, alt: "n")[$n$] in the following way: #math.equation(block: true, alt: "Y sub n open parenthesis x close parenthesis equals lim sub p → n the fraction cos open parenthesis p π close parenthesis J sub p open parenthesis x close parenthesis minus J sub minus p open parenthesis x close parenthesis over sin open parenthesis p π close parenthesis .")[$Y_(n) ( x ) = lim_(p → n) frac(cos ( p π ) J_(p) ( x ) − J_(− p) ( x ), sin ( p π )) .$] As each linear combination of #math.equation(block: false, alt: "J sub p open parenthesis x close parenthesis")[$J_(p) ( x )$] and #math.equation(block: false, alt: "J sub minus p open parenthesis x close parenthesis")[$J_(− p) ( x )$] is a solution to Bessel's equation of order #math.equation(block: false, alt: "p")[$p$], then as we take the limit as #math.equation(block: false, alt: "p")[$p$] goes to #math.equation(block: false, alt: "n")[$n$], #math.equation(block: false, alt: "Y sub n open parenthesis x close parenthesis")[$Y_(n) ( x )$] is a solution to Bessel's equation of order #math.equation(block: false, alt: "n")[$n$]. It also turns out that #math.equation(block: false, alt: "Y sub n open parenthesis x close parenthesis")[$Y_(n) ( x )$] and #math.equation(block: false, alt: "J sub n open parenthesis x close parenthesis")[$J_(n) ( x )$] are linearly independent. Therefore when #math.equation(block: false, alt: "n")[$n$] is an integer, we have the general solution to Bessel's equation of order #math.equation(block: false, alt: "n")[$n$] #math.equation(block: true, alt: "y equals A J sub n open parenthesis x close parenthesis plus B Y sub n open parenthesis x close parenthesis ,")[$y = A J_(n) ( x ) + B Y_(n) ( x ) ,$] for arbitrary constants #math.equation(block: false, alt: "A")[$A$] and #math.equation(block: false, alt: "B")[$B$]. Note that #math.equation(block: false, alt: "Y sub n open parenthesis x close parenthesis")[$Y_(n) ( x )$] goes to negative infinity at #math.equation(block: false, alt: "x equals 0")[$x = 0$]. Many mathematical software packages have these functions #math.equation(block: false, alt: "J sub n open parenthesis x close parenthesis")[$J_(n) ( x )$] and #math.equation(block: false, alt: "Y sub n open parenthesis x close parenthesis")[$Y_(n) ( x )$] defined, so they can be used just like say #math.equation(block: false, alt: "sin open parenthesis x close parenthesis")[$sin ( x )$] and #math.equation(block: false, alt: "cos open parenthesis x close parenthesis")[$cos ( x )$]. In fact, they have some similar properties. For example, #math.equation(block: false, alt: "minus J sub 1 open parenthesis x close parenthesis")[$− J_(1) ( x )$] is a derivative of #math.equation(block: false, alt: "J sub 0 open parenthesis x close parenthesis")[$J_(0) ( x )$], and in general the derivative of #math.equation(block: false, alt: "J sub n open parenthesis x close parenthesis")[$J_(n) ( x )$] can be written as a linear combination of #math.equation(block: false, alt: "J sub n minus 1 open parenthesis x close parenthesis")[$J_(n − 1) ( x )$] and #math.equation(block: false, alt: "J sub n plus 1 open parenthesis x close parenthesis")[$J_(n + 1) ( x )$]. Furthermore, these functions oscillate, although they are not periodic. See Figure #math.equation(block: false, alt: "1")[$1$] for graphs of Bessel functions. #figure(figph[Two side-by-side plots for x from 0 to 10, each a pair of decaying oscillations about a quarter period out of step. On the left the blue curve starts at 1 and the green at 0, peaking near 0.58; on the right both plunge toward minus infinity as x approaches 0.], alt: "Two side-by-side plots for x from 0 to 10, each a pair of decaying oscillations about a quarter period out of step. On the left the blue curve starts at 1 and the green at 0, peaking near 0.58; on the right both plunge toward minus infinity as x approaches 0.", caption: [Figure #math.equation(block: false, alt: "1")[$1$]: Plot of the #math.equation(block: false, alt: "J sub 0 open parenthesis x close parenthesis")[$J_(0) ( x )$] and #math.equation(block: false, alt: "J sub 1 open parenthesis x close parenthesis")[$J_(1) ( x )$] in the first graph and #math.equation(block: false, alt: "Y sub 0 open parenthesis x close parenthesis")[$Y_(0) ( x )$] and #math.equation(block: false, alt: "Y sub 1 open parenthesis x close parenthesis")[$Y_(1) ( x )$] in the second graph.]) #examplebox("Example 4")[Using Bessel functions to Solve a ODE][ Other equations can sometimes be solved in terms of the Bessel functions. For example, given a positive constant #math.equation(block: false, alt: "λ")[$λ$], #math.equation(block: true, alt: "x y double prime plus y prime plus λ squared x y equals 0 ,")[$x y^(″) + y^(′) + λ^(2) x y = 0 ,$] can be changed to #math.equation(block: false, alt: "x squared y double prime plus x y prime plus λ squared x squared y equals 0")[$x^(2) y^(″) + x y^(′) + λ^(2) x^(2) y = 0$]. Then changing variables #math.equation(block: false, alt: "t equals λ x")[$t = λ x$] we obtain via chain rule the equation in #math.equation(block: false, alt: "y")[$y$] and #math.equation(block: false, alt: "t")[$t$]: #math.equation(block: true, alt: "t squared y double prime plus t y prime plus t squared y equals 0 ,")[$t^(2) y^(″) + t y^(′) + t^(2) y = 0 ,$] which can be recognized as Bessel's equation of order 0. Therefore the general solution is #math.equation(block: false, alt: "y open parenthesis t close parenthesis equals A J sub 0 open parenthesis t close parenthesis plus B Y sub 0 open parenthesis t close parenthesis")[$y ( t ) = A J_(0) ( t ) + B Y_(0) ( t )$], or in terms of #math.equation(block: false, alt: "x")[$x$]: #math.equation(block: true, alt: "y equals A J sub 0 open parenthesis λ x close parenthesis plus B Y sub 0 open parenthesis λ x close parenthesis .")[$y = A J_(0) ( λ x ) + B Y_(0) ( λ x ) .$] This equation comes up for example when finding fundamental modes of vibration of a circular drum, but we digress. ] === Footnotes \[1\] Named after the German mathematician Ferdinand Georg Frobenius (1849 – 1917). \[2\] See Joseph L. Neuringera, The Frobenius method for complex roots of the indicial equation, International Journal of Mathematical Education in Science and Technology, Volume 9, Issue 1, 1978, 71–77. \[3\] Named after the German astronomer and mathematician Friedrich Wilhelm Bessel (1784 – 1846).