# Message Boards

Answer
(Unmark)

GROUPS:

Piecewise series solutions to ODEs Frobenius and power series solutions combine accuracy and speed. However, they generally only have a finite convergence radius, which limits their applications. A solution to this issue is to “sew” several series solutions together into a piecewise function. Piecewise solutions retain series solutions’ accuracy and speed while having a potentially unlimited domain. The current essay starts with presenting computational tools for finding the analytical formulas for Frobenius and power series solutions to linear homogeneous ODEs; these tools can be used to compute approximate power and Frobenius series solutions. After this, a function to compute piecewise series solutions is presented. It automatically picks the parameters such that the error of the solution it outputs doesn’t exceed some fixed tolerance. The essay then lists and addresses the possible issues. The essay is a successor to another post on the same topic. The algorithm for finding Frobenius series solutions given in the present essay is about 15 times faster than the old one; while the old algorithm does substitution explicitly, like a human, the new algorithm uses a prepared formula. The current essay also presents a better algorithm for finding piecewise solutions - it is way faster and requires less input from the user. In[]:= SaveDefinitions True;
Quick demonstration of the results Compute the formula for the Frobenius series solution to the ODE 2 2 2 x x=0 frobeniusDSolveFormula[2 2 x 2 (x+3) Out[]= - 1 2 1 3 (-12+10n-2 2 n 2 r 2 n 2 r -3+21(n+r)+18(-1+n+r)(n+r) The first part output tells us that r -1/2 1/3 a n y[x]= ∞ ∑ n=0 a n r+n x y[x]≈ N ∑ n=0 a n r+n x N=40 In[]:= exampleApproxFrobSolns=frobeniusDSolve[2 2 x 2 (x+3) Plot[#,{x,0,2},ImageSizeMedium,AxesLabel{"x","y[x]"}]&/@exampleApproxFrobSolns Out[]= , Because the ODE has another singular point at x=-3 x=3 In[]:= exampleApproxPwSolns=frobeniusPiecewiseDSolve[2 2 x 2 (x+3) Legended[Row[{Plot[{exampleApproxPwSolns[[1]],exampleApproxFrobSolns[[1]]},{x,0,6},ImageSizeMedium,PlotRange{-4,4},AxesLabel{"x","y[x]"}]," ",Plot[{exampleApproxPwSolns[[2]],exampleApproxFrobSolns[[2]]},{x,0,6},ImageSizeMedium,PlotRange{0,3},AxesLabel{"x","y[x]"}]}],LineLegend[{ColorData[97][1],ColorData[97][2]},{"Piecewise","Frobenius"}]] Out[]=
Piecewise series solutions are very accurate while also being fast. Below is given an approximate solution to 2 2 2 x x=0 2 2 2 x -12 10 69 Out[]=
The functions presented in the essay can also be used to solve ODEs of higher order, or ODEs that don’t have a singular point. Below, for example, are the four linearly independent solutions to (4) x In[]:= powerSolns=frobeniusDSolve[x''''[t]+x[t]0,x,t,0,40]; Plot[#,{t,-6,6},ImageSizeSmall]&/@powerSolns Out[]= , , ,
Introduction Given a linear homogenous ODE (ordinary differential equation) with polynomial coefficients, we can always write it in the form M ∑ m=0 Z ∑ p=0 C mp m (x- x 0 p d p dx ( 1 ) where x 0 C x= x 0 u[t]= ∞ ∑ n=0 a n r+n t ( 2 ) where t≡x- x 0 0= ∞ ∑ k=-Z r+k t k+Z ∑ n=Max[0,k-M] L kn a n ( 3 ) where the formula for L kn t k+Z ∑ n=Max[0,k-M] L kn a n ( 4 ) This gives us s a set of equations that allows us to determine the constant r a n C L kn Min[Z,M+n-k] ∑ p=Max[0,n-k] C k-n+p,p n Π i=n+1-p ( 5 )
Code to extract the matrix C To obtain a Frobenius solution, we want to extract the matrix C t In[]:= ClearAll[changeVarsToT];Protect[FrobTvar];changeVarsToT[eqn_, depvar_, indvar_, x0_] := eqn /. depvar[indvar] depvar[FrobTvar], Derivative[p_][depvar][indvar] Derivative[p][depvar][FrobTvar] /. { indvar FrobTvar + x0 }; changeVarsToT[y''[x] 2 x Out[]= ′′ y 2 FrobTvar We then transform an ODE of the form y''[t]-y[x]-2y'[x] {y''[x],y[x],2y'[x]} y''[x]+y[x]+2y'[x]0 In[]:= ClearAll[extractTerms];extractTerms[eqn_Equal] := Apply[List, Expand@eqn[[1]] - Expand@eqn[[2]]]; extractTerms[ ′′ y Out[]= {y[t],2 ′ y ′′ y Each term has the form m t p d p dt ,m,p In[]:= ClearAll[containsQ]; (* test if expr contains any of the symbols syms *)containsQ[expr_, syms_] := If[AtomQ[expr], AnyTrue[syms, SameQ[expr, #]&], containsQ[Head[expr], syms] || AnyTrue[Apply[List, expr], containsQ[#, syms]&] ]; In[]:= ClearAll[mpform];mpform[term_, depvar_Symbol] := Block[ factors = If[SameQ[Head[term], Times], Apply[List, term], {term}], = 1, m = 0, p = Null , Which[ MatchQ[#, Derivative[p1_][depvar][FrobTvar]], If[SameQ[p, Null], p = # /. {Derivative[p1_][depvar][FrobTvar] p1}, p = $Failed], MatchQ[#, depvar[FrobTvar]], If[SameQ[p, Null], p = 0, p = $Failed], MatchQ[#, m1_ FrobTvar m1_ FrobTvar mpform[#,y]&/@{y[FrobTvar],2 ′ y ′′ y Out[]= {{1,0,0},{2,0,1},{1,0,2}} The output tells us that there are 3 terms, and their numbers {,m,p} {1,0,0} {2,0,1} {1,0,2} y[t]+2y'[t]+y''[t]=0 C In[]:= ClearAll[extractCmatrix];extractCmatrix[eqn_Equal, depvar_, indvar_, x0_] := Block[ terms, mplist, Cmatrix, M, Z , terms = extractTerms@changeVarsToT[eqn, depvar, indvar, x0]; mplist = mpform[#, depvar]& /@ terms; If[MemberQ[mplist, $Failed], {$Failed, $Failed, $Failed}, M = Max[mplist[[All, 2]]]; Z = Max[mplist[[All, 3]]]; Cmatrix = ConstantArray[0, {1 + M, 1 + Z}]; (Cmatrix[[1 + #[[2]], 1 + #[[3]]]] += #[[1]])& /@ mplist; {M, Z, Cmatrix} ]]; Here is what the function outputs for the equation 3xy''[x]+y[x]0 M,P,C M,P MP C extractCmatrix[3xy''[x]+y[x]0,y,x,0] Out[]= {1,2,{{1,0,0},{0,0,3}}} MatrixForm[{{1,0,0},{0,0,3}}] Out[]//MatrixForm=
The matrix has two nonzero entries. The first one, C 00 y[x] C 12 3xy''[x]
Formula for the Frobenius solution
The recursion relation For M≤k k+Z ∑ n=k-M L kn a n ( 6 ) Given that L kn n=k+Z k-M+R ∑ n=k-M L kn a n where R≤Z+M L kn n=k-M+R a n a n 1 L n+M-R,n n-1 ∑ k=n-R L n+M-R,k a k ( 7 ) We don’t have to worry about L n+M-R,n R R L k+M-R,k L k+M-R,k Min[Z,R] ∑ p=Max[0,R-M] C M-R+p,p k Π i=k+1-p In order for this expression to be zero independently of r C R C M-R+p,p p Max[0,R-M] Min[Z,R] R R=Z+M C R R-1 C R R R wrong L k+M-R,k In[]:= ClearAll[computeLsmatrix];computeLsmatrix[M_, Z_, Cmatrix_] := Table[ Sum[ Cmatrix[[1 + k - n + p, 1 + p]] Product[r + i, {i, n + 1 - p, n}], {p, Max[0, n - k], Min[Z, M + n - k]} ], {k, -Z, M - 1}, {n, 0, Z + M - 1}] In[]:= ClearAll[findRecursionReln];findRecursionReln[M_, Z_, Cmatrix_, L_] := Block { R }, R = Z + M; While[AllTrue[(Cmatrix[[1 + M - R + #, 1 + #]] 0)& /@ Range[Max[0, R - M], Min[Z, R]], #&], R--; ]; R, a[n] 1 L[n + M - R, n]
The indicial equation For -Z≤k<M k+Z ∑ n=0 L kn a n ( 8 ) If we use the fact that L kn n>k Z+M a (Z+M)(Z+M) L s L s L M=3 Z=2 Out[]//MatrixForm=
When using the method of Frobenius, the indicial equation is generally determined by looking at the coefficient before the lowest power of t k -Z L -Z,0 a 0 a 0 L -Z,0 C 0,Z 0 Π i=1-Z ( 9 ) If C 0,Z (r-(Z-1))(r-(Z-2))...r=0 r∈{0,1,2,...Z-1} t u[t]= ∞ ∑ k=0 a k r+k t u C 0,Z C 1,Z 2 t C 2,Z |