Message Boards

GROUPS:

Staff Picks Mathematics Physics Calculus Equation Solving Wolfram Language Numerical Computation Symbolic Computations Computer-Based Maths

Rauan Kaldybaev

Piecewise series solutions to ODEs

Rauan Kaldybaev

Posted 9 months ago

1670 Views

1 Reply

4 Total Likes

Follow this post

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

(x+3)

y''+7x(x+3)y'-3y=0

near the singular point

x=0

frobeniusDSolveFormula[2

(x+3)

y''[x]+7x(x+3)y'[x]-3y[x]0,y,x,0]

Out[]=

-

,{(-3+21(1+r)+18r(1+r))a[1]-7r-12(-1+r)r,(-3+21(2+r)+18(1+r)(2+r))a[2]-2(-1+r)r-(7(1+r)+12r(1+r))a[1],(-3+21(3+r)+18(2+r)(3+r))a[3]-2r(1+r)a[1]-(7(2+r)+12(1+r)(2+r))a[2]},a[n]

(-12+10n-2

+10r-4nr-2

)a[-2+n]+(-17+29n-12

+29r-24nr-12

)a[-1+n]

-3+21(n+r)+18(-1+n+r)(n+r)



The first part output tells us that

can be either

-1/2

1/3

, and the rest allows us to determine the coefficients

y[x]=

∞

∑

n=0

r+n

. Compute the approximate solution with a finite number of terms,

y[x]≈

∑

n=0

r+n

, with

N=40

In[]:=

exampleApproxFrobSolns=frobeniusDSolve[2

(x+3)

y''[x]+7x(x+3)y'[x]-3y[x]0,y,x,0,40];

Plot[#,{x,0,2},ImageSizeMedium,AxesLabel{"x","y[x]"}]&/@exampleApproxFrobSolns

Out[]=





Because the ODE has another singular point at

x=-3

, the Frobenius series solutions start diverging at

x=3

. The piecewise solution, however, goes on being accurate:

In[]:=

exampleApproxPwSolns=frobeniusPiecewiseDSolve[2

(x+3)

y''[x]+7x(x+3)y'[x]-3y[x]0,y,x,0,{0,6}];

Legended[Row[{Plot[{exampleApproxPwSolns[[1]],exampleApproxFrobSolns[[1]]},{x,0,6},ImageSizeMedium,PlotRange{-4,4},AxesLabel{"x","y[x]"}]," ",Plot[{exampleApproxPwSolns[[2]],exampleApproxFrobSolns[[2]]},{x,0,6},ImageSizeMedium,PlotRange{0,3},AxesLabel{"x","y[x]"}]}],LineLegend[{ColorData[97][1],ColorData[97][2]},{"Piecewise","Frobenius"}]]

Out[]=

	Piecewise
	Frobenius

Piecewise series solutions are very accurate while also being fast. Below is given an approximate solution to

(x+3)

y''+7x(x+3)y'-3y=0

near

x=0

. Its error, measured as

2

(x+3)

y''+7x(x+3)y'-3y

, does not exceed

-12

. The solution only took

milliseconds to compute.

Out[]=

Approximate solution to -3y[x]+7x(3+x)

′

[x]+2

(3+x)

′′

[x] = 0

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)

[t]+x[t]=0

In[]:=

powerSolns=frobeniusDSolve[x''''[t]+x[t]0,x,t,0,40];

Plot[#,{t,-6,6},ImageSizeSmall]&/@powerSolns

Out[]=





Introduction

Given a linear homogenous ODE (ordinary differential equation) with polynomial coefficients, we can always write it in the form

∑

m=0

∑

p=0

(x-

)

(

)

where

is some (any) number and

is a matrix with constant entries. To solve this equation using the method of Frobenius near

, assume solution of the form

u[t]=

∞

∑

n=0

r+n

(

)

where

t≡x-

. Substituting this into the equation (1) yields

∞

∑

k=-Z

r+k

k+Z

∑

n=Max[0,k-M]

(

)

where the formula for

is given below. To satisfy the equation, the coefficients before all powers of

must be zero:

k+Z

∑

n=Max[0,k-M]

=0-Z≤k

(

)

This gives us s a set of equations that allows us to determine the constant

and coefficients

in (2). In a sense, the matrix

encapsulates all the relevant information about the ODE (1) that we need to compute a solution.

≡

Min[Z,M+n-k]

∑

p=Max[0,n-k]

k-n+p,p

i=n+1-p

(r+i)

(

)

Code to extract the matrix
C

To obtain a Frobenius solution, we want to extract the matrix

. For that, we first change the independent variable form whatever it used to be to

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]

y[x],y,x,0]

Out[]=

′′

[FrobTvar]

FrobTvar

y[FrobTvar]

We then transform an ODE of the form

y''[t]-y[x]-2y'[x]

into a list of terms

{y''[x],y[x],2y'[x]}

. The meaning of this list is that

y''[x]+y[x]+2y'[x]0

In[]:=

ClearAll[extractTerms];extractTerms[eqn_Equal] := Apply[List, Expand@eqn[[1]] - Expand@eqn[[2]]];

extractTerms[

′′

[t]-y[t]-2y'[x]]

Out[]=

{y[t],2

′

[x],

′′

[t]}

Each term has the form



; let us compute the numbers

,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

], m = # /.

m1_

FrobTvar

 m1, MatchQ[#, FrobTvar], m = 1, True, If[SameQ[Head[#], depvar] || containsQ[#, {FrobTvar, depvar}],  = $Failed,  *= #] ]& /@ factors; If[IntersectingQ[{Null, $Failed}, {m, p}], $Failed, {, m, p} ]];

mpform[#,y]&/@{y[FrobTvar],2

′

[FrobTvar],

′′

[FrobTvar]}

Out[]=

{{1,0,0},{2,0,1},{1,0,2}}

The output tells us that there are 3 terms, and their numbers

{,m,p}

are

{1,0,0}

{2,0,1}

, and

{1,0,2}

, respectively, which corresponds to the ODE

y[t]+2y'[t]+y''[t]=0

. Now, we can write a function to compute the matrix

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

. The output is of the form

M,P,C

, where

M,P

are the dimensions of the rectangular

MP

matrix

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=



1	0	0
0	0	3



The matrix has two nonzero entries. The first one,

, corresponds to the term

y[x]

, and the other one,

, corresponds to

3xy''[x]

Formula for the Frobenius solution

The recursion relation

For

M≤k

, the equation (4) becomes

k+Z

∑

n=k-M

(

)

Given that

would often be zero for

n=k+Z

or other high values, we can further simplify the above equation by writing

k-M+R

∑

n=k-M

where

R≤Z+M

is the highest number under the condition that

is nonzero for the upper limit of summation

n=k-M+R

. From here, we obtain the recursion relation that tells us

in terms of the previous coefficients:

n+M-R,n

n-1

∑

k=n-R

n+M-R,k

(

)

We don’t have to worry about

n+M-R,n

being identically zero anymore. Now, we need to ask the question of how to find

. The number

is defined by the condition that

k+M-R,k

is nonzero. The formula (5) tells us that

k+M-R,k

Min[Z,R]

∑

p=Max[0,R-M]

M-R+p,p

i=k+1-p

(r+i)

In order for this expression to be zero independently of

, all constants

must be zero. Then,

is the highest number such that at least one of the coefficients

M-R+p,p

is nonzero, where

goes from

Max[0,R-M]

Min[Z,R]

. One way to compute

is the following. Step 1: start with

R=Z+M

. Step 2: check if the

’s are all zero. Step 3: If true change

R-1

and return to step 2; otherwise, which is if any of the

’s is nonzero, stop. Because we start from the highest possible value of

and go down in steps of 1, the algorithm is going to return the correct

and not some other

wrong

that passes the criterion

k+M-R,k

≠0

but isn’t the highest number to do so.

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] 

L[n + M - R, n]

Sum[Expand[-L[n + M - R, k]]a[k], {k, n - R, n - 1}];

The indicial equation

For

-Z≤k<M

, the equation (4) can be written as

k+Z

∑

n=0

=0-Z≤k<M

(

)

If we use the fact that

for

n>k

, we can read (8) as a matrix equation. The

Z+M

dimensional vector

lies in the nullspace of the

(Z+M)(Z+M)

matrix

, where

is a “slice” of the larger matrix

. Below you can see the first few entries of the matrix for

M=3

Z=2

. Notice that it’s in lower-triangular form.

Out[]//MatrixForm=

(-1+r)r C 02	0	0	0	0
r C 01 +(-1+r)r C 12	r(1+r) C 02	0	0	0
C 00 +r C 11 +(-1+r)r C 22	(1+r) C 01 +r(1+r) C 12	(1+r)(2+r) C 02	0	0
C 10 +r C 21 +(-1+r)r C 32	C 00 +(1+r) C 11 +r(1+r) C 22	(2+r) C 01 +(1+r)(2+r) C 12	(2+r)(3+r) C 02	0
C 20 +r C 31	C 10 +(1+r) C 21 +r(1+r) C 32	C 00 +(2+r) C 11 +(1+r)(2+r) C 22	(3+r) C 01 +(2+r)(3+r) C 12	(3+r)(4+r) C 02

When using the method of Frobenius, the indicial equation is generally determined by looking at the coefficient before the lowest power of

. For the lowest value of

, which is

-Z

, we have

-Z,0

. Since we demand

to be nonzero, we have

-Z,0

. Substituting (4) yields

0,Z

i=1-Z

(r+i)=0

(

)

0,Z

≠0

, then

(r-(Z-1))(r-(Z-2))...r=0

, and thus

r∈{0,1,2,...Z-1}

. This means that the powers of

in the formula

u[t]=

∞

∑

k=0

r+k

(2) are all integers, and we effectively get a power series solution. This makes sense because the coefficient before the highest derivative of

0,Z

1,Z

2,Z

+...

, and

Piecewise series solutions to ODEs

Quick demonstration of the results

Introduction

Code to extract the matrix C

Formula for the Frobenius solution

The recursion relation

The indicial equation

Code to extract the matrix
C