The equation can be put in invariant reduced form as a biconfluent type

Heun equation

y''-(x^2+1-k/x+1/(4x^2))y=0.

Then it is possible to apply Kovacic algorithm to this equation.

In fact, the solvability by Liouvillian functions was studied by A. Duval and M.Loday-

Richaud using Kovacic's algorithm(see proposition 13 of the attached file).

If you are seeking a numerical method, you should first determine the type of ODE, initial or boundary value problem, and then specify the initial/boundary conditions exactly. In this way there are some possibilities to solve the problem numerically.

But if you are trying to find the exact solution, I suggest the Frobenius series solutions which can find in any differential equation textbook.

You can find a solution using Maple software. The answer is then in terms of Heun functions and their integrals.

The equation can be put in invariant reduced form as a biconfluent type

Heun equation

y''-(x^2+1-k/x+1/(4x^2))y=0.

Then it is possible to apply Kovacic algorithm to this equation.

In fact, the solvability by Liouvillian functions was studied by A. Duval and M.Loday-

Richaud using Kovacic's algorithm(see proposition 13 of the attached file).

We have published a Matlab tool box just for this task. It is available at:

http://www.mathworks.at/matlabcentral/fileexchange/41354

The theory behing this code can be found at:

http://arxiv.org/abs/1304.3312

We are presently preparing a significant update to this tool box to improve the results when solving Sturm-Liouville problems

regards

Paul

If you don't get anywhere, please contact me, I believe if you give me the values for the constrants, the range in which the solution is required and the boundary conditions to make the solution unique then I can deliver numerical solutions without must work.

I concur with Davoud Mirzaei in regards to the additional info you need to pose this problem well. I am not sure how you obtained this equation for Ficks law (the flux has the direction of the negative gradient), but assuming that x is the independent variable and y the dependent variable, then the following trick makes this problem almost trivial: let z = y'. Hence y'' = z'. Then you end up with a system of uncoupled ODEs:

(1) z' = -(2x^2+1) *z/x,

(2) y'= z.

Solve (1) to obtain z= exp(-x^2)/x. Plug this into (2) and you obtain y = Ei(-x^2)/2 + C, where Ei is the Exponential Integral function. This is a fast answer, so please excuse any typos, mistakes, and/or omissions. Also, I am not sure what your variables mean, so what I just wrote might not apply to this problem (you did not provide enough information). In any case, it is interesting.

I would like to thank everyone for their very helpful suggestions, especially Angelo and Paul for their suggestion to use Heun integrals. I would also like to provide an additional background into the problem:

I wish to solve Fick's Law (cylindrical system) on an infinite domain to obtain the flux on a ring (with k-fold symmetry). I begin with my governing equations:

dC/dt = D[d2C/dr^2+ (1/r)dC/dr + (1/r^2)d2C/dtheta^2 ]

where D is the diffusivity ( a constant). The boundary conditions are:

C=1 at t=0 for all r and theta

C=0 at r=1 for all t>0, for all theta

C=1 at r--> \infty (infinity)

dC/dtheta =0 at theta=2*pi/k (where k is an integer)

I assume the form of the solution to be given by:

C= C'(r,t)*cos(k*theta)

which leads to:

dC'/dt=D[d2C'/dr^2+ (1/r)dC'/dr - (k/r)^2C' ]

Using a similarity transformation eta = r/sqt(4Dt), I obtain the following ODE:

d2C'/deta^2 +(2*eta+1/eta)*dC'/deta-(k/eta)^2*C'=0

This is the form of a second-order ODE with variable coefficients. If possible, I would wish to avoid obtaining the solution as an infinite series (Frobenius). I am unsure how to formulate the problem as a Heun equation. My problem does not seem to be of the typical Heun or Lame type, is it?

MAPLE that's enough.

Dear Ricardo, Maple is NOT enough. While a CAS is helpful, mathematicians cannot simply become technicians operating a sophisticated instrument. *Somebody* has to program Maple. Besides, Maple is unable to find solutions in a very large number of cases.

Is it not possible to solve this equation using series method? since the coefficients have only low order powers of x, the recurrence equations in the series method should not be very complicated.

You can transform the differential equation to a linear differential system of dimension 2

(Moser reduction algorithm) and then use the ISOLDE on MAPLE to get the solution.

Have you tried power series solution?

Divyaraj Desai

If it is enough to solve numerically, I can saggest one old method. You can obtaine Rikkati equations for u=y'/y. It is very easy and stable to solve it and no problen with boundary value problem. Also You can esealy abtaine asimptotic for big x.

Only You must controle the direction of numerically solving to take into account right exponent (the same way as in WKB).

Good Luck

Sorry for some missprint

You can find series solution of above diffrentail eq. x=0 is a regular singular point therefore we can use Frobenious method to find solution.

Hi once again. If I am correct you wish to solve the equation:

d2C'/deta^2 +(2*eta+1/eta)*dC'/deta-(k/eta)^2*C'=0

Do you have initial values or boundary values for C and a value for k? If yes then there is a simple numerical solution to this problem. If you can send me the initial conditions and k I will do the computation and deliver the matlab code for the solution.

I can deliver a solution as a power series with respect to the Gram polynomials.

regards Paul

See the Trefethen's book spectral method in matlab. You will find every thing about ODE and PDE in this book.

Here is the link.

Best regards,

http://www.mathworks.ir/downloads/english/Spectral%20Methods%20in%20MATLAB.pdf

Utiliser Runge-Kutta si c'est un problème à conditions initiales (EDO du second ordre: 2 conditions initiales y(0) et sa dérivée en 0).

Si le problème est un problème aux conditions aux limites, utiliser la méthode du tir (Shooting method; Reference: Ha, S.N.: A nonlinear shooting method for two-point boundary value problems. Computers and Mathematics with Applications 42(10-11), 1411-1420 (2001)) + Runge-Kutta.

K=0, la solution est une équation intégrale

y(x)=c1+c2*Intégrale[Exp[-2*u+1/u; u variant de 1 à x]

c1 et c2 des constantes.

We may use power series method to solve this type of ODE. More specifically,

Frobenius method may be applied to find the solution.

We can use popular numerical method of Runga Kutta 4th order method for multi ODE or multidimentional pde.

There are numerous methods to solve the second order linear equation what have already been discussed above but you use some computing software like matlab th open source version is scilab

With proper substitution , the second order ODE can be simplified. Next step can be followed as per the characteristics of the problem. for example finding eigen value a better solution can be obtained.

To solve 2nd order ODE, apply series solution method.

i.e. solution of bessel and legendre differential equations these methods are available in most of all text/reference books of differential equations

E. Kreyszig: Advanced Engineering Mathematics, 8th Edition,Jhon Wiley & Sons

It would be convenient to obtain an analytical solution (closed-form or series solution)to an ODE with variable coefficients but a series solution can be limited by its convergence domain. Every ODE can be solved numerically provided the convergence and stability criteria are established. If the problem is an initial value one whereby y(a) and y'(a) are prescribed then the best solution method is the Runge Kutta 4th order method or predictor-corrector method (Euler-Cauchy or Heun). If the problem is a boundary value problem with y(a) and y(b) prescribed then the ODE should be tackled using finite difference method (FDM) OR Finite element method (FEM). The later method is more versatile, converges faster and is to be preferred. However, the FEM or ODE's is not understood by a substantial number of researchers. Later I will provide the FE formulation of the ODE using the elementary piecewise linear shape function [x/l 1-x/l] ; l = b-a. This simple shape function can be surprisingly accurate!

The link below is spreadsheet to solve a steady state one dimensional ground water flow problem in a vertical column. The governing equation is a 2nd order ODE with variable coefficients solved by finite difference and both matrix solver and iterative built in Excel solver. Hope it may help.

https://www.researchgate.net/publication/235955596_A_spreadsheet_model_to_solve_steady_state_groundwater_flow_equation_in_a_non-homogeneous_soil_column

Another model that simulate diffusion by Ficks law is random walk theory, below is another spreadsheet simulates particles diffusion in a box by random walk method.

https://www.researchgate.net/publication/258032344_Simulation_of_Brownian_Motion_on_ExcelSheet

Data A spreadsheet model to solve steady state groundwater flow e...

Data Simulation of Brownian Motion on ExcelSheet

Such second order IVP/BVP can easily be solved by applying one of the methods:

1) The method of the superposition

2) The method of chasing

3) The Adjoint operator methods

4) Finite difference approximation procedure

Surprisingly Divyaraj Desai has not stated whether he is solving an IVP or BVP!

Dear Amaechi Anyaegbunam,

It is don't mension!!!

Let us consider Lv=0 where L second order linear operator. With general condition

av(x_0)+bv'(x_0)=1 and cv(l)+ev'(l)=0 any condition can be reduced to that one. More!

Let us change v(x)=u(x)/(1-a1u(x_0)_b1u(x_0) then we have the same equation for u(x) (it depend from v by multiplication on const) and new condition

(a+a1)u(x_0)+(b+b1)u'(x_0)=1. Now we can choose a+a1=1 and b+b1=0.

Now we have

Lu=0, u(x_0)=1, au(l)+bu'(l)=0.

Ok! Let us introduse w(x)=u'(x)/u(x). We have Ricatti eqution for w with initial condition at l so w(l)=-a/b. We can solve IVP for Ricatti eqution from l to x_0. As we know w(x) we have u=exp(int from x_0 to x w(y) dy).

You can see that u(x_0)=1 and u'(l) =u(l)w(l)exp(...) so

au(l)+bu(l)(-a/b)=0 it is identity!

Such a way You can solve boundary value problen by solving 2 initialy valied problem of first order.

In the case of eigen value problem the way more conplex but simple enough.

In the case of Green function no problem also.

Good Luck!

Some notions about Ricatti eqution. We have no problem with singularity in coefficients or solution, we have no problem with numerical stability and accuracy. More it is not needed differentiation of coefficents. You have stochatis ode with white noise in coefficints? No problem!

My dear Constantin Koshel

It matters a lot whether it is an IVP or BVP if the problem is to be tackled numerically. See my previous comments.

Also, it matters if a problem is an IVP or BVP when a solution is sought for a PDE like the heat equation.

Some people have suggested that a series solution can be obtained for the given ODE using Frosbenius method. They seem to forget that x = 0 is an irregular singular point for the ODE such that Frosbenius method fails because x.p(x) is infinite at x = 0 where p(x) = (2*x^2+1)/x^2.

Dear Constantin. You have mentioned one method out of many methods of solving the ODE. Your preference may not be to my liking. Could you please offer to us the analtytic solution to the given ODE: x^2.y'' +(2x^2+1)y' -ky =0

Dear Amaechi Anyaegbunam,

1. I say that numerically BVP and IVP are simple both (not a diference) for second order linear ODE (not PDE). There are many methods to establish that fact, I indicate one.

2. Analitical theory of linear ODE give us a way to solve ODE as series (with log maltiplaer) no problem with singularite. It give us generaal solution!!! may be conversation not so good.

3. The ricatti equation have no problem with singularity when solving numerically).

4. Only very resrtrict case of second order linear ODE have analitical solution (see Kamke), in elemantary functions. The solution as spetial function can be found (may be hypergeometric :))

Good Luck.

About xp(x) is infinite. You can solve Euler equation

x^ny^(n)+ax^n-1y^(n-1) + ... ets.?

The solution contains all singularities You need. :)

Good Luck

My dear Constantin please don't be too rash in rushing out answers. The 2nd order Euler equation is ao*x^2*y"(x) + a1*x*y'(x) + a2*y(x) = 0. In this case p(x) = a1/(ao*x) so that x. p(x) = a1/ao is finite at x = 0. Thus, the 2nd order Euler equation has a regular singular point at x = 0 !

I am trying to solve variable coefficients ODE by converting into constant coefficients

Amaechi Anyaegbunam

Sorry,

I do not understand what You meen.

aox^2y''+a1xy'+a2y=0

we can put a1=0?

so we have

y''+(a2/ao)(1/x^2)y

the same case as x^2y''+(2x^2+1)y'-ky=0.

here we have (1/x^2)q(x) where q(x)=sum x^i

the first aproximation we can take from Euler equation

y''+a(1/x^2)y=0

y=C1*x^g1+C2*x^g2 where gi are solutions of

g(g-1)+a=0.

It can be singular if g

Dear Amaechi,

WRT the IVP or BVP query, the unsteady state problem is defined as a IVBVP, which was reduced to a BVP using a similarity transformation [x=r/sqt(4t)].

Dear All,

You can see the following link for getting the arthematic equation from the differential equation

http://www.stewartcalculus.com/data/CALCULUS%20Concepts%20and%20Contexts/upfiles/3c3-2ndOrderLinearEqns_Stu.pdf

http://www.math.psu.edu/tseng/class/Math251/Notes-2nd%20order%20ODE%20pt1.pdf

for solving using any numerical methods,

http://people.maths.ox.ac.uk/suli/nsodes.pdf

http://faculty.olin.edu/bstorey/Notes/DiffEq.pdf

http://isc.temple.edu/physics402/notes/iv_notes.pdf

Best Regards,

Ravikiran

Dear Constantin I am worried that we are arguing over an elementary criterion for solvability in series expansion of a 2nd order ODE. The 2nd order ODE y"(x)+p(x)y'(x) +q(x)y(x) is said to have a regular singular point at x=0 if both x*p(x) and x^2*q(x) are both finite at x=0, where '*" stands for multiplication. If any of these criteria are not satisfied then x=0 is an irregular singular point and the Frosbenius method cannot be used to obtain a series solution. Thus the Euler equation ao*x^2*y"(x) + a2*y(x) = 0 has a regular singular point at x = 0. This is my last comment on this issue!

Please, may I clarify myself. When I say the solution of an ODE cannot be obtained in series form I do not mean it has no closed-form solution or that its solution cannot be obtained numerically.

Dear Amaechi Anyaegbunam,

Sorry, I not Good understand You fist comment, so the initial eqution have the regular singular poit

(2x^2+1)/x^2=p(x) and x^2p(x)=1+2x^2 and x^2p(x) =1 when x=0, so You first comment not connected with that equation, I do not understand it.

Sorry again!

Good Luck!

Konstantin Koshel

Second order linear homogenous ODE is in form of Cauchy-Euler S form or Legender form you can convert it in to linear with constant coefficient ODE which can solve by standard methods.But variable coefficients are not as above than most general method is call vatiation of parameter methods by using lemma you can find two independent solution and with help of wronskian you can find general solution.If highest derivative is multiplied by small parameter and singularities aries than use method of match asymptomatic expansions or composite expansions method given in book of "singular perturbation techniques".by ALI Nayfeh.wielly publication,1974. or Omali .I hope all 2nd ODE can solve from above hints.Use any std books of differential equations and their solutionAsk for any kind of Nonliner PDE BVP or IVP of diffussion type..

qestion of Amaechi, singular perturbation technique will be useful.Above books may be useful.

I send my answer as a file attachment.doc. Your sincerely Anna tomova

Dear Manoj Mehta,

May be it is good idea to give here some best books list (Nayef as example). When I try to find student book for ODE where the equations f(u',u,x) was considered (for case cann't be represented in form u'=g(u,x)) I can find only 1967 Year book.

Good Luck

Ali Nayfeh is for sinular perturbation.But Differential equtions and their solutions byMurphy gives direct solution of any kind of ODE and PDE.Method of variation of parameter is most general method to find solutions of linear ode with variable coefficient.

Dear Divyaraj Desai

y(x)= ∑_(n=0)^∞▒a_(n ) 〖(x-α)〗^(n ) α≠0

y^(' ) (x)=∑_(n=1)^∞▒〖na_(n ) 〖(x-α)〗^(n-1) 〗

y^('' ) (x)=∑_(n=2)^∞▒〖n(n-1)〖(x-α)〗^(n-2) 〗

x^(2 ) y^('' ) (x)+(2x^(2 )+1) y^(' ) (x)-ky(x)=0

a_(n )=?(a_(0 ) ,a_1 ) n≥2

Dear Nasir Taghizadeh,

I would make a probe to answer Your question for the formulas of the coefficients, but we aren't sure if the series converge or not. For Riccati's equation, if the functions of the variable x are absolute convergent (as in the this case), there exists always a solution in series which are convergent . Please, answer me.Your sincerely, Anna Tomova.

Dear Anna Tomova,

We can calculare the radius of convegence by usual way. Then we can use some transforms and obtain solutions as series in another interval. But it seems that You are right, the methods of Riccati eqution more simle and convinient here.

Good Luck!

I am not happy with errors occurring in answers provided by research gate members. In the paper by Duval and Laday-Richaud provided by Juan Morales-Ruiz it is shown that using the transformation y=u.exp(Integral(bdx/2/a) the ODE a(x)y"+b(x)y'+c(x)y=0 can be converted to the reduced form u'' -r(x)u = 0 where

r(x) = (2ab'-2ba'+b^2-4ac)/4/a^2.

Could Morales-Ruiz explain to us how the ODE x^2y" +(2x^2+1)y'-ky = 0 transforms into the ODE y" - (x^2+1-k/x+1/4/x^2)y = 0 when we substitute a = x^2, b = 2x^2+1, c= -k in the expression for r(x)?

Dear Amaechi Anyaegbunam,

You are right!

Good Luck

Here attached a matlab code to solve this ODE. Best regards

with cheb.m from Trefethen book : Spectral methods in

matlab. Best regards.

% Solve ODE f1(x)*u_xx + f2(x)*u_x + (x)*u = f4(x), with

%Dirichlet or/and Neumann boundary conditions

clear all

close all

N = 20;

% differentiation matrix (see trefethen's book : Spectral Method in

% Matlab

[D,x] = cheb(N);

D2 = D^2;

D1=D;

I=eye(N+1);

O=0*I;

% find the boundaries

bl = find(x==-1);

br = find(x==1);

% fonctions of the ODE

f1 = x.*x;

f2 = 2.*x.*x+1;

f3=3.;

f4=0.*x;

% ODE operator

L=diag(f1)*D2+diag(f2)*D1+diag(f3)*I;

% Boundary conditions

% Dirichlet for x=-1

L(bl,:) = O(bl,:); L(bl,:) = I(bl,:);

% Neumann for x=-1

%L(bl,:) = O(bl,:); L(bl,:) = D1(bl,:);

%Dirichlet for x= +1

%L(br,:) = O(br,:); L(br,:) = I(br,:);

%Neumann for x=+1

L(br,:) = O(br,:); L(br,:) = D1(br,:);

f4(bl)=-20; % u(-1)=-20

f4(br)=5; % u_x(+1)=5

%solution

u = L\f4;

%plot the solution

clf, subplot('position',[.1 .4 .8 .5])

plot(x,u,'.','markersize',16)

xx = -1:.01:1;

uu = polyval(polyfit(x,u,N),xx);

line(xx,uu)

xlabel x, ylabel y

Dear Slim Kaddeche. Thanks for the code you have provided. I will now test its accuracy with the Bessel ODE x^2y" + xy' +(x^2-k^2)y = 0. If it is correct you have provided the first definite answer to Divyaraj Desai's question provided of course he is solving a BVP with either/or Dirichlet and Neumann boundary conditions.

I had argued about the solvability in series expansion of the ODE x^2y" +(2x^2+1)-ky = 0. I hereby attach extracts from the textbook by Stephenson to support my opinion that no series expansion can be obtained. Apologies for the poor quality of the extract from the old 1973 textbook.

Manoj Mehta mentioned that the method of variation of parameters is a general method for finding solutions of linear ODE with variable coefficients. This statement is not strictly true. The method of variation of parameters is used for finding the particular solutions (particular integrals) of the non-homogeneous ODE a(x)y"+b(x)y'+c(x)y = f(x) provided you have succeeded in obtaining one solution of the homogeneous problem a(x)y"+b(x)y'+c(x)y = 0. Divyaraj Desai is looking for the solution to a homogeneous ODE and the method of variation of parameters cannot help.

Assuming what you need is the numerical solution...

Personally, I would start by rewriting this in state-space form:

x^2*y''+(2x^2+1)*y'-ky=0 ->

f(x)*y'' + g(x)*y' - k*y = 0

State-Space solution:

y1 = y

y2 = y1' = y'

y2' = y'' = -g(x)/f(x)*y2 + k/f(x)*y1

Then, the state space matrix equation is y'(x) = (y1'(x) y2'(x))T = A(x)*y(x)

where

y1' = y2

y2' = -g(x)/f(x)*y2 + k/f(x)*y1

Therefore, A(x) =

0 1

k/f(x) -g(x)/f(x)

Then your next step is to solve this fun nonlinear equation...

If you only need the solution for small displacements about a point, you can linearize about that point and your life will be easy:

http://mmae.iit.edu/~mpeet/Classes/MMAE443/443Lecture03.pdf

If you need to solve this globally, you need to consider other methods to solve it as-is in nonlinear form:

http://en.wikibooks.org/wiki/Control_Systems/Time_Variant_System_Solutions

Consult a control systems textbook if you're still having issues! Good luck!

Dear Divyaraj Desai

y(x)= ∑_(n=0)^∞▒a_(n ) 〖(x-α)〗^(n ) α≠0

y^(' ) (x)=∑_(n=1)^∞▒〖na_(n ) 〖(x-α)〗^(n-1) 〗

y^('' ) (x)=∑_(n=2)^∞▒〖n(n-1)〖(x-α)〗^(n-2) 〗

x^(2 ) y^('' ) (x)+(2x^(2 )+1) y^(' ) (x)-ky(x)=0

a_(n )=?(a_(0 ) ,a_1 ) n≥2

Dear Nasir Taghizadeh,

Сolleague Amaechi Anyaegbunam some times have argued that initial equation

x^2y''+(2x^2+1)y'-ky=0.

have no regular solution at x=0!

We do not hear it!

The equation y''-(x^2+1-k/x+1/(4x^2))y=0. not equovalent first one!

I suggest simple furmula to solve the first equation by recurent procedure.

y=uexp(fi)

u(i+1)=u(i)*(Hi*cosh(Hi*h)+vi*sinh(Hi*h))/Hi

and

v(i+1)=Hi*(Hi*sinh(Hi*h)+vi*cosh(Hi*h))/(Hi*cosh(Hi*h)+vi*sinh(Hi*h))

where h- step of integration (can be hi) and Hi=f(h*i+h/2)

and f(x)=sqrt(-(1/x^3-1/(4x^4))-1-(1+k)/x^2)

see attachment for details.

Good Luck

I see it is not fastest way but very stable.

It seems yhat singularite in solution at list exp(5/4x) at x->0

Thу algorithm in general form for f(x) and good for example for Bessel ODE :)

Dear Colleague Divyaraj Desai,

Decompose the coefficients and solve them as a first order ODE as in the theorem. The strategy also works for the higher order ODE.

Dear Colleague Divyaraj Desai

You can solve it via the Homotopy perturbation method which give you an approximated solution or simply solve it using Mathematica package using the command Dsolve.

Dear colleague and Mummin Mohammed Hussein do you know that the Homotopy perturbation method can fail to give a correct solution when applied to ODEs. See the article by Francisco M. Fernandez arXiv:0808.2078

You can solve it by Power series method. One can also reduce the variable coefficient problem into constant coefficient problem by a suitable transformation or one can solve by using numerical methods or simple software's like MATLAB or by MATHEMATICA. For example you can look at the follwoing links

1. https://math.stackexchange.com/questions/646558/second-order-homogeneous-differential-equation-with-non-constant-coefficients

2. http://www.math.ualberta.ca/~xinweiyu/334.1.10f/DE_2nd_order_eq.pdf

page no.22