# Message Boards

Answer
(Unmark)

GROUPS:

Frobenius and power series solutions to ODEs in Wolfram Not all ordinary differential equations, or ODEs, can be solved analytically, and power series or Frobenius series solutions can then be extremely helpful. The current essay provides Wolfram Language function to help automate power and Frobenius series solutions, a manual work that is done much better by computer than by human. The essay looks at linear homogenous ODEs with polynomial coefficients. A major disadvantage of series solutions is that generally, they only converge within a limited region, which severely restricts their usage; in this essay, a method for piecewise Frobenius/power series is given, that allows to overcome the challenge of limited convergence radius by sewing Frobenius series solution with several power series solutions. In[]:= SaveDefinitions True;
Quick demonstration of the results Solve the equation 2 2 x x=0 In[]:= frobeniusDSolve[2 2 x
Out[]= -3, 1 2 7(-1+n+r)a[-1+n] -3+7(n+r)+2(-1+n+r)(n+r) The output tells us that r -3 1/2 a n 7(-1+n+r) a n-1 -3+7(n+r)+2(-1+n+r)(n+r) In[]:= frobeniusDSolve[f''[x]+f[x],f,x,0] Out[]= {0,1},{r(1+r)a[1]0},a[n]- a[-2+n] (-1+n+r)(n+r) r 0 x 1 x 2 2 x 10 In[]:= frobeniusNDSolve[2 2 x Out[]= 1- 21x 5 49 2 x 5 343 3 x 15 3 x x 7x 18 49 2 x 264 1715 3 x 20592 16807 4 x 494208 117649 5 x 9335040 823543 6 x 193489920 117649 7 x 89303040 823543 8 x 2190901248 5764801 9 x 57994444800 40353607 10 x 1648263168000 In[]:= {Plot[%[[1]],{x,0,2},ImageSizeMedium],Plot[%[[2]],{x,0,2},ImageSizeMedium]} Out[]= , The more terms we compute, the more accurate the solution becomes. Here is what the second solution ( x 7x 18 49 2 x 264 Out[]=
With ncoefs = 100, it is possible to reduce the error to be of order -15 10 Out[]=
While the ‘regular’ Frobenius solution seems to start diverging at x=3 x=-3 2 2 x 2 (x+3) Out[]=
Code used to generate the figures
Step-by-step solutions To automate Frobenius’s method, let us solve an ODE manually and automate each step of the process as we go. Consider the following ODE: df dx ( 1 ) It has no singular points, but for illustrative purposes that’s fine. The first step is to assume solution of the form f[x]= r (x- x 0 ∞ ∑ n=0 a n n (x- x 0 ( 2 ) where x 0 t≡x- x 0 u[t]≡f[x- x 0 t≡x- x 0 x 0 ( 3 ) The equation (1) has no singular points, so, just for the sake of example, let us look for a series solution near x=1 x 0 t u u'+tu=0 ( 4 ) where u'[t]≡ du dt In[]:= ClearAll[frobeniusVarRename];frobeniusVarRename[eqn_Plus, depvar_, indvar_, singularPoint_] := Block[ {t, u, p}, Expand[eqn /. depvar[indvar] u[t], Derivative[p_][depvar][indvar] Derivative[p][u][t] /. indvar t + singularPoint ] ] In[]:= frobeniusVarRename[f'[x]+(x-1)f[x],f,x,1] Out[]= tu[t]+ ′ u Let us proceed to the second step. In terms of t u u[t]= ∞ ∑ n=0 a n r+n t ( 5 ) Plugging it into (4), we get ∞ ∑ n=0 a n r+n-1 t ∞ ∑ n=0 a n r+n t ( 6 ) From here, we can determine the coefficients a n In[]:= ClearAll[frobeniusSubstitute];frobeniusSubstitute[eqn_] := Block[ {a, r, p, q, n, u, t}, eqn /. u[t] Sum[a[n] r + n t q_ t r + n + q t r + n - p t r + n + 1 - p t q_ t r + n + q - p t ∞ ∑ n=0 a n r+n t In[]:= frobeniusSubstitute[tu[t]+ ′ u Out[]= ∞ ∑ n=0 -1+n+r t ∞ ∑ n=0 n+r t To compute the coefficients a ∑ n r+n x t ∞ ∑ n=-1 n+r t ∞ ∑ n=1 n+r t ( 7 ) To perform this step, the function below first decomposes the left-hand side of the equation into a list of terms. For each term of the form ∞ ∑ n=0 a n r+n+ n 0 t a n r+n+ n 0 t n 0 ∞ ∑ n= n 0 a n- n 0 r+n t ∞ ∑ n= n 0 a n- n 0 r+n t In[]:= ClearAll[frobeniusNRSumForm];frobeniusNRSumForm[eqn_] := Block terms = eqn[[#]]& /@ Range[Length[eqn]], (* transforms a + b + c 0 into {a, b, c} *) summand, n0, p, r, n, t , Table summand = term /. {Sum Function[{expr, limits}, expr]}; (* finds the summand in ∞ ∑ n=0 a n r + n + n 0 t p_ t p_ t n 0 ∞ ∑ n=0 a n r + n + n 0 t a n r + n + n0 t ∞ ∑ n= n 0 a n-n0 r + n t ∞ ∑ n=0 -1+n+r t a n ∞ ∑ n=0 n+r t a n ∞ ∑ n=-1 n+r t a n+1 ∞ ∑ n=1 n+r t a n-1 In[]:= frobeniusNRSumForm ∞ ∑ n=0 -1+n+r t ∞ ∑ n=0 n+r t Out[]= ∞ ∑ n=-1 n+r t ∞ ∑ n=1 n+r t The indicial equation for the equation is obtained by demanding the coefficient before the lowest power of t n min r a 0 ( 8 ) From here, r=0 n=1 n=0 (1+r) a 1 ( 9 ) The “bottom equations”, which are the equations (7) for small n n≥1 t n 0 ∞ ∑ n= n 0 n r+n t n 0 ( n 0 min ( n 0 max In[]:= ClearAll[frobeniusBottomEqns];frobeniusBottomEqns[nrsform_] := Block[ r, n, t, n0list, n0min, n0max , n0list = #[[2, 2]]& /@ nrsform; n0min = Min[n0list]; n0max = Max[n0list]; Table[ Coefficient[ Total[If[#[[2, 2]] ≤ currn, #[[1]] /. {n currn}, 0 ]& /@ nrsform], r + currn t ∞ ∑ n=-1 n+r t a n+1 ∞ ∑ n=1 n+r t a n-1 a 0 a 1 In[]:= frobeniusBottomEqns ∞ ∑ n=-1 n+r t ∞ ∑ n=1 n+r t Out[]= {ra[0]0,(1+r)a[1]0} For all higher values of n ∑ n a n+1 a n-1 r+n t t (1+n+r) a n+1 a n-1 Substituting the solution r=0 a n a n-2 n+r ( 10 ) For n≥2 In[]:= ClearAll[frobeniusRecursionReln];frobeniusRecursionReln[nrsform_] := Block t, n, r, a, terms, factors, n0max, newterms, oldcoefs, newcoefs, getfactors = Function[term, If[SameQ[Head[term], Times], Table[ term[[i]], {i, Length[term]} ], term ]] , terms = Coefficient[nrsform[[#, 1]], n + r t a n a n a n-1 Total[oldcoefs] Coefficient[Total[newcoefs], a[n]] In[]:= frobeniusRecursionReln ∞ ∑ n=-1 n+r t ∞ ∑ n=1 n+r t Out[]= a[n]- a[-2+n] n+r We can now define a function that would put together what we’ve done so far. In[]:= ClearAll[frobeniusDSolve, frobeniusDSolveStepByStepHelper];Options[frobeniusDSolve] = {"StepByStep" False};frobeniusDSolve[eqn_, depvar_, indvar_, singularPoint_, OptionsPattern[]] := Block[ frobform, nrsform, bottomeqns, recre |