2.1.782 Problem 804

Solved as second order ode using Kovacic algorithm
Maple
Mathematica
Sympy

Internal problem ID [9954]
Book : Collection of Kovacic problems
Section : section 1
Problem number : 804
Date solved : Wednesday, March 05, 2025 at 08:01:23 AM
CAS classification : [[_Emden, _Fowler]]

Solve

xy+3y+4x3y=0

Solved as second order ode using Kovacic algorithm

Time used: 0.330 (sec)

Writing the ode as

(1)xy+3y+4x3y=0(2)Ay+By+Cy=0

Comparing (1) and (2) shows that

A=x(3)B=3C=4x3

Applying the Liouville transformation on the dependent variable gives

z(x)=yeB2Adx

Then (2) becomes

(4)z(x)=rz(x)

Where r is given by

(5)r=st=2AB2BA+B24AC4A2

Substituting the values of A,B,C from (3) in the above and simplifying gives

(6)r=16x4+34x2

Comparing the above to (5) shows that

s=16x4+3t=4x2

Therefore eq. (4) becomes

(7)z(x)=(16x4+34x2)z(x)

Equation (7) is now solved. After finding z(x) then y is found using the inverse transformation

y=z(x)eB2Adx

The first step is to determine the case of Kovacic algorithm this ode belongs to. There are 3 cases depending on the order of poles of r and the order of r at . The following table summarizes these cases.

Case

Allowed pole order for r

Allowed value for O()

1

{0,1,2,4,6,8,}

{,6,4,2,0,2,3,4,5,6,}

2

Need to have at least one pole that is either order 2 or odd order greater than 2. Any other pole order is allowed as long as the above condition is satisfied. Hence the following set of pole orders are all allowed. {1,2},{1,3},{2},{3},{3,4},{1,2,5}.

no condition

3

{1,2}

{2,3,4,5,6,7,}

Table 2.782: Necessary conditions for each Kovacic case

The order of r at is the degree of t minus the degree of s. Therefore

O()=deg(t)deg(s)=24=2

The poles of r in eq. (7) and the order of each pole are determined by solving for the roots of t=4x2. There is a pole at x=0 of order 2. Since there is no odd order pole larger than 2 and the order at is 2 then the necessary conditions for case one are met. Since there is a pole of order 2 then necessary conditions for case two are met. Therefore

L=[1,2]

Attempting to find a solution using case n=1.

Looking at poles of order 2. The partial fractions decomposition of r is

r=4x2+34x2

For the pole at x=0 let b be the coefficient of 1x2 in the partial fractions decomposition of r given above. Therefore b=34. Hence

[r]c=0αc+=12+1+4b=32αc=121+4b=12

Since the order of r at is Or()=2 then

v=Or()2=22=1

[r] is the sum of terms involving xi for 0iv in the Laurent series for r at . Therefore

[r]=i=0vaixi(8)=i=01aixi

Let a be the coefficient of xv=x1 in the above sum. The Laurent series of r at is

(9)r2ix3i16x39i1024x727i32768x11405i4194304x151701i134217728x1915309i8589934592x2372171i274877906944x27+

Comparing Eq. (9) with Eq. (8) shows that

a=2i

From Eq. (9) the sum up to v=1 gives

[r]=i=01aixi(10)=2ix

Now we need to find b, where b be the coefficient of xv1=x0=1 in r minus the coefficient of same term but in ([r])2 where [r] was found above in Eq (10). Hence

([r])2=4x2

This shows that the coefficient of 1 in the above is 0. Now we need to find the coefficient of 1 in r. How this is done depends on if v=0 or not. Since v=1 which is not zero, then starting r=st, we do long division and write this in the form

r=Q+Rt

Where Q is the quotient and R is the remainder. Then the coefficient of 1 in r will be the coefficient this term in the quotient. Doing long division gives

r=st=16x4+34x2=Q+R4x2=(4x2)+(34x2)=4x2+34x2

We see that the coefficient of the term x in the quotient is 0. Now b can be found.

b=(0)(0)=0

Hence

[r]=2ixα+=12(bav)=12(02i1)=12α=12(bav)=12(02i1)=12

The following table summarizes the findings so far for poles and for the order of r at where r is

r=16x4+34x2

pole c location pole order [r]c αc+ αc
0 2 0 32 12

Order of r at [r] α+ α
2 2ix 12 12

Now that the all [r]c and its associated αc± have been determined for all the poles in the set Γ and [r] and its associated α± have also been found, the next step is to determine possible non negative integer d from these using

d=αs()cΓαcs(c)

Where s(c) is either + or and s() is the sign of α±. This is done by trial over all set of families s=(s(c))cΓ until such d is found to work in finding candidate ω. Trying α=12 then

d=α(αc1)=12(12)=0

Since d an integer and d0 then it can be used to find ω using

ω=cΓ(s(c)[r]c+αcs(c)xc)+s()[r]

The above gives

ω=(()[r]c1+αc1xc1)+()[r]=12x+()(2ix)=12x2ix=12x2ix

Now that ω is determined, the next step is find a corresponding minimal polynomial p(x) of degree d=0 to solve the ode. The polynomial p(x) needs to satisfy the equation

(1A)p+2ωp+(ω+ω2r)p=0

Let

(2A)p(x)=1

Substituting the above in eq. (1A) gives

(0)+2(12x2ix)(0)+((12x22i)+(12x2ix)2(16x4+34x2))=00=0

The equation is satisfied since both sides are zero. Therefore the first solution to the ode z=rz is

z1(x)=peωdx=e(12x2ix)dx=eix2x

The first solution to the original ode in y is found from

y1=z1e12BAdx=z1e123xdx=z1e3ln(x)2=z1(1x3/2)

Which simplifies to

y1=eix2x2

The second solution y2 to the original ode is found using reduction of order

y2=y1eBAdxy12dx

Substituting gives

y2=y1e3xdx(y1)2dx=y1e3ln(x)(y1)2dx=y1(ie2ix24)

Therefore the solution is

y=c1y1+c2y2=c1(eix2x2)+c2(eix2x2(ie2ix24))

Will add steps showing solving for IC soon.

Maple. Time used: 0.006 (sec). Leaf size: 21
ode:=diff(diff(y(x),x),x)*x+3*diff(y(x),x)+4*x^3*y(x) = 0; 
dsolve(ode,y(x), singsol=all);
 
y=c1sin(x2)+c2cos(x2)x2

Maple trace

`Methods for second order ODEs: 
--- Trying classification methods --- 
trying a quadrature 
checking if the LODE has constant coefficients 
checking if the LODE is of Euler type 
trying a symmetry of the form [xi=0, eta=F(x)] 
checking if the LODE is missing y 
-> Trying a Liouvillian solution using Kovacics algorithm 
   A Liouvillian solution exists 
   Group is reducible or imprimitive 
<- Kovacics algorithm successful`
 

Maple step by step

Let’s solve(d2dx2y(x))x+3ddxy(x)+4x3y(x)=0Highest derivative means the order of the ODE is2d2dx2y(x)Isolate 2nd derivatived2dx2y(x)=4x2y(x)3(ddxy(x))xGroup terms withy(x)on the lhs of the ODE and the rest on the rhs of the ODE; ODE is lineard2dx2y(x)+3(ddxy(x))x+4x2y(x)=0Check to see ifx0=0is a regular singular pointDefine functions[P2(x)=3x,P3(x)=4x2]xP2(x)is analytic atx=0(xP2(x))|x=0=3x2P3(x)is analytic atx=0(x2P3(x))|x=0=0x=0is a regular singular pointCheck to see ifx0=0is a regular singular pointx0=0Multiply by denominators(d2dx2y(x))x+3ddxy(x)+4x3y(x)=0Assume series solution fory(x)y(x)=k=0akxk+rRewrite ODE with series expansionsConvertx3y(x)to series expansionx3y(x)=k=0akxk+r+3Shift index usingk>k3x3y(x)=k=3ak3xk+rConvertddxy(x)to series expansionddxy(x)=k=0ak(k+r)xk+r1Shift index usingk>k+1ddxy(x)=k=1ak+1(k+1+r)xk+rConvertx(d2dx2y(x))to series expansionx(d2dx2y(x))=k=0ak(k+r)(k+r1)xk+r1Shift index usingk>k+1x(d2dx2y(x))=k=1ak+1(k+1+r)(k+r)xk+rRewrite ODE with series expansionsa0r(2+r)x1+r+a1(1+r)(3+r)xr+a2(2+r)(4+r)x1+r+a3(3+r)(5+r)x2+r+(k=3(ak+1(k+1+r)(k+r+3)+4ak3)xk+r)=0a0cannot be 0 by assumption, giving the indicial equationr(2+r)=0Values of r that satisfy the indicial equationr{2,0}The coefficients of each power ofxmust be 0[a1(1+r)(3+r)=0,a2(2+r)(4+r)=0,a3(3+r)(5+r)=0]Solve for the dependent coefficient(s){a1=0,a2=0,a3=0}Each term in the series must be 0, giving the recursion relationak+1(k+1+r)(k+r+3)+4ak3=0Shift index usingk>k+3ak+4(k+4+r)(k+6+r)+4ak=0Recursion relation that defines series solution to ODEak+4=4ak(k+4+r)(k+6+r)Recursion relation forr=2ak+4=4ak(k+2)(k+4)Solution forr=2[y(x)=k=0akxk2,ak+4=4ak(k+2)(k+4),a1=0,a2=0,a3=0]Recursion relation forr=0ak+4=4ak(k+4)(k+6)Solution forr=0[y(x)=k=0akxk,ak+4=4ak(k+4)(k+6),a1=0,a2=0,a3=0]Combine solutions and rename parameters[y(x)=(k=0akxk2)+(k=0bkxk),a4+k=4ak(k+2)(4+k),a1=0,a2=0,a3=0,b4+k=4bk(4+k)(k+6),b1=0,b2=0,b3=0]
Mathematica. Time used: 0.052 (sec). Leaf size: 41
ode=x*D[y[x],{x,2}]+3*D[y[x],x]+4*x^3*y[x]==0; 
ic={}; 
DSolve[{ode,ic},y[x],x,IncludeSingularSolutions->True]
 
y(x)4c1eix2ic2eix24x2
Sympy. Time used: 0.196 (sec). Leaf size: 20
from sympy import * 
x = symbols("x") 
y = Function("y") 
ode = Eq(4*x**3*y(x) + x*Derivative(y(x), (x, 2)) + 3*Derivative(y(x), x),0) 
ics = {} 
dsolve(ode,func=y(x),ics=ics)
 
y(x)=C1J12(x2)+C2Y12(x2)x