0
$\begingroup$

I want to solve below ODE:

$\partial_x y(x) = \frac{a \: g(x)}{5 x^2}+ \frac{c}{x^2}$

In which $a$ and $c$ are constants and $g(x)$ is as below:

g[x_] := 1/(11664 a^4)*E^(-(((5 + 3 Sqrt[30]) a)/(35 x))) (36 (-5 c E^(((5 + 3 Sqrt[30]) a)/(35 x)) (324 a^2 - 630 a x - 7595 x^2) + 324 a^4 (E^((6 Sqrt[6/5] a)/(7 x)) (18 a^2 - 21 Sqrt[30] a x + 245 x^2) C1 + (18 a^2 + 21 Sqrt[30] a x + 245 x^2) C2)) + 245 c (E^((6 Sqrt[6/5] a)/(7 x)) (18 Sqrt[30] a^2 - 630 a x + 245 Sqrt[30] x^2) ExpIntegralEi[((5 - 3 Sqrt[30]) a)/(35 x)] - (18 Sqrt[30] a^2 + 630 a x + 245 Sqrt[30] x^2) ExpIntegralEi[((5 + 3 Sqrt[30]) a)/(35 x)]));

C1 = -((2790 c + 5832 a^4 C2 - 245 Sqrt[30] c ArcCoth[3Sqrt[6/5]])/(5832 a^4));

C2 =-(5 c (44100 (434 + 67 a) E^(1/35 (5 + 3 Sqrt[30]) a) + 2 E^(6/7 Sqrt[6/5]a) (558 (-17150 + a (245 (-5 + 6 Sqrt[30]) + 3 a (35 (-18 + Sqrt[30]) + 6 (-5 + 3 Sqrt[30]) a))) + 245 (3430 Sqrt[30] +  a (245 (-36 + Sqrt[30]) + 18 a (-35 + 21 Sqrt[30] - 18 a + Sqrt[30] a))) ArcCoth[3 Sqrt[6/5]]) + 245 (3430 Sqrt[30] + a (245 (-36 + Sqrt[30]) + 18 a (7 (-5 + 3 Sqrt[30]) + (-18 + Sqrt[30]) a))) E^( 6/7 Sqrt[6/5] a) ExpIntegralEi[1/35 (5 - 3 Sqrt[30]) a] - 245 (3430 Sqrt[30] + a (245 (36 + Sqrt[30]) + 18 a (35 + 21 Sqrt[30] + 18 a +  Sqrt[30] a))) ExpIntegralEi[1/35 (5 + 3 Sqrt[30]) a]))/(11664 a^4 (17150 + 1225 a + 1470 Sqrt[30] a + 1890 a^2 + 105 Sqrt[30] a^2 + 90 a^3 + 54 Sqrt[30]a^3 + (-17150 + a (245 (-5 + 6 Sqrt[30]) + 3 a (35 (-18 + Sqrt[30]) + 6 (-5 + 3 Sqrt[30]) a))) E^(6/7 Sqrt[6/5] a)));

I used below code, but it can not be run:

Simplify[DSolve[yy'[x] == a a20/5 g[x]/x^2 + c/x^2, yy[x], x ]]

Could anyone help me to solve the equation? The boundary condition for $y$ are:

$\partial_xy(x=1)=-s$

and

$y(x=\infty)=w$

which $w$ and $s$ are known constants.

If there is no analytic solution: $c$ is a constant which should be find by a boundary condition. Also, values for $a$ should be chosen larger than 1 in magnitude. For example $a=2$ or $a=-2$.

$\endgroup$
8
  • $\begingroup$ psitwo is not defined. $\endgroup$
    – anderstood
    Feb 22, 2018 at 18:56
  • $\begingroup$ Sorry, it was a typo. I edit it @anderstood $\endgroup$
    – Wat Watson
    Feb 22, 2018 at 18:57
  • $\begingroup$ DSolve returns {{yy[x] -> -(a20/(58320 a^3 x)) - c/x + C[1]}}. C1 and C2 are not used. $\endgroup$
    – anderstood
    Feb 22, 2018 at 19:02
  • $\begingroup$ C1 and C2 are ysed in definition of g(x). and this solution is not correct because $g$ is a function of $x$ @anderstood $\endgroup$
    – Wat Watson
    Feb 22, 2018 at 19:05
  • $\begingroup$ OK I see, I edited your code so that the copy and paste now gives the right expression for g. Would you be OK with a numerical solution instead of a closed-form (which might not exist)? That would require values for a,b,c. $\endgroup$
    – anderstood
    Feb 22, 2018 at 19:11

1 Answer 1

0
$\begingroup$

Assuming s=w=a20=1 the following works:

c = c /. First@Solve[(a*g[x]/(5 x^2) + c/x^2 /. x -> 1 // N) == -s, c]
sol[x_] = 
NDSolveValue[yy'[x] == a/5*g[x]/x^2 + c/x^2 && yy[10^10] == w, yy[x], {x, 1, 10}]

Plot[sol[x], {x, 1, 10}]

You can update with your actual values that you forgot to indicate.

Regarding the closed-form solution: I think you ODE is too complicated for a closed-form to exist, or to be found by Mathematica. But that's just a belief, not a proof...

$\endgroup$
4
  • $\begingroup$ I put s = 1; w = 100; a = 2; and I plot the solution in LogLogPlot but y-axis of the output is not logarithmic. Do you know where is the problem? LogLogPlot[sol[x], {x, 1, 100}] $\endgroup$
    – Wat Watson
    Feb 22, 2018 at 19:56
  • $\begingroup$ If you have a new question, feel free to start a new post. (Try LogLogPlot[sol[x], {x, 1, 100}, PlotRange -> Full]) $\endgroup$
    – anderstood
    Feb 22, 2018 at 20:00
  • $\begingroup$ I think there is a problem in your answer. When I use below code, boundary condition at infinity (here 500) is not satisfied in the output: s = -1; w = 1000; a = 2; c = c /. First@ Solve[(a*g[x]/(5 x^2) + c/x^2 /. x -> 1 // N) == -s, c] sol[x_] = NDSolveValue[yy'[x] == a/5*g[x]/x^2 + c/x^2 && yy[500] == w, yy[x], {x, 1, 10}] Plot[sol[x] , {x, 1, 500}, PlotRange -> Full] $\endgroup$
    – Wat Watson
    Feb 22, 2018 at 20:29
  • $\begingroup$ @WatWatson In NDSolveValue, you need to adjust {x, 1, 10} to the range of interest (here, {x, 1, 500}). $\endgroup$
    – anderstood
    Feb 22, 2018 at 20:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.