7
$\begingroup$

The ODE I need to solve is

$$\left(y^3y^{\prime\prime\prime}\right)^\prime+\frac{5}{8}xy^\prime-\frac{1}{2}y+\frac{c}{y}=0$$ where $\prime$ denotes differentiation, $c$ is a constant and $0<c\le1$. It needs 4 boundary conditions (BCs): $y^\prime(0)=y^{\prime\prime\prime}(0)=0$ are given, and the other two values of $y(0)$ and $y^{\prime\prime}(0)$ are determined by shooting for BCs at infinity (or $x_\text{Max}$).

Inspired by @ bbgodfrey's idea in this answer, I know that the far-field asymptotic behavior is given by DSolveing the remainder of the ODE after the nonlinear derivative term are neglected:

DSolve[5/8 x y'[x] - 1/2 y[x] + c/y[x] == 0, y[x], x]

Thus the far-field asymptotic behavior is $y= kx^{4/5}$, where $k$ is a constant. Accordingly, we refer to a solution with $y\sim x^{4/5}$ asymptotic behavior as a $x^{4/5}$ solution. Here, I want to search for this kind of solutions by shooting method in which the ODE is integrated from $x_0=0$ to some large value of $x_\text{Max}$, say $x_\text{Max}=10$.

The shooting parameters $y(0)$ and $y^{\prime\prime}(0)$ will be adjusted so that $y-\frac{5}{4}xy^\prime=0$ or $y+\frac{25}{4}x^2y^{\prime\prime}=0$ at $x=x_\text{Max}$ is satisfied. The two BCs were chosen to be independent to involve low-order derivatives and to require $y\propto x^{4/5}$ at the end of the integration interval.

Thanks to @ xzczd's answer to this question, an example of the code is:

c = 1;
(*c=1/10;*)
ode = D[y[x]^3*y'''[x], x] + 5/8 x y'[x] - 1/2 y[x] + c/y[x] == 0;
bc1 = y[x] - 5/4 x y'[x] == 0;
bc2 = y[x] + 25/4 x^2 y''[x] == 0;

Thank again to this answer, I can obtain the series solution at $x_0=0$ with seriesDSolve

solSeries = seriesDSolve[ode, y, {x, 0, 5}, 
                        {y[0] -> a, y'[0] -> 0, y''[0] -> b, y'''[0] -> 0}]

We thus obtain a set of boundary conditions at $x=x_0$ by

newbclist =  Thread[(Derivative[#][y][x0] == 
             (D[Normal@solSeries, {x, #}] /. x -> x0) &) /@ Range[0, 3]]

Then, give the interval, solve the ODE and plot the parameter plane

x0 = 0; xMax = 6;
sol = ParametricNDSolveValue[{ode, newbclist}, y, {x, x0, xMax}, {a, b}]
ContourPlot[{bc1, bc2} /. x -> xMax /. y -> sol[a, b] // Evaluate, 
{a, 0.2, 0.6}, {b, -0.05, 1}]

Unfortunately, Mathematica gives something like this:

ParametricNDSolveValue::ndsz: At x$803007 == 0.27144285125721357`, step size is effectively zero; singularity or stiff system suspected.

as encounters in this problem.

After some experiments, I found the situation appears to be worse for a larger parameter $c$, i.e. $c\rightarrow 1$. And I also played with the intervals of $a$ and $b$ but it seems no help. Three ugly results for your reference (sorry!...)

1. When I use c = 1/10, xMax = 0.5 with {a, -0.2, 1}, {b, -0.1, 1}, ContourPlot gives

enter image description here

2. When I use c = 1/10, xMax = 0.5 with {a, 0.2, 0.6}, {b, -0.05, 1}, ContourPlot gives

enter image description here

3. When I change to bc1 = 1/2 y[x] - 5/8 x y'[x] - c/y[x] == 0; with c = 1, xMax = 6 and {a, -0.2, 1}, {b, -0.1, 1}, ContourPlot gives

enter image description here

To sum up, my questions are:

(1) After I find the parameter $y(0)$ and $y^{\prime\prime}(0)$, how can I figure out the associated constant $k$ in the far-field asymptotic behavior $y= kx^{4/5}$?

(2) Is it possible to estimate the parameters $y(0)$ and $y^{\prime\prime}(0)$ automatically to start the integration, as shown in this answer?

Improve this question
$\endgroup$
2
  • $\begingroup$ I think it's just the nature of your equation: Some choices of {a, b} can't lead to a solution extending to infinity. Just observe this: mid = sol[0.6, 0.6]; {{lb, rb}} = mid["Domain"]; Plot[mid[t], {t, lb, rb}, PlotRange -> All] $\endgroup$
    – xzczd
    Commented Jan 16, 2016 at 4:29
  • $\begingroup$ Thanks, @ xzczd, That really is the question. I am wondering is there any technique to find out the ranges of $a$ and $b$ so that the solutions can shoot up to an appropriate $x_\text{Max}$. In other words, I just try to explore the existence of some isolated parameter pairs of $(a,b)$ for certain given $c$ and $x_\text{Max}$. That is, the values of $y(0)$ and $y''(0)$ are determined by shooting for boundary condition at $x_\text{Max}$. $\endgroup$
    – W. Robin
    Commented Jan 16, 2016 at 6:56

1 Answer 1

Reset to default
2
$\begingroup$

Automatically looking for proper $y(0)$ and $y''(0)$ isn't easy, at least I can't think out a solution, still, I don't think you need such fascinating technique to solve the new equation, you just need to slightly modify the code for plotting to filter out those improper $(a,b)$:

Quiet@ContourPlot[
  If[y["Domain"][[1, -1]] < xMax, 1, #] == 0 & /@ Subtract @@@ {bc1, bc2} /. 
     x -> xMax /. y -> sol[a, b] // Evaluate, {a, 1, 1.6}, {b, -0.1, 0.2}, 
  PlotPoints -> 50]

enter image description here

The plot range is obtained through not too much trial and error.

"How can I figure out the associated constant $k$ in the far-field asymptotic behavior $y=kx^{4/5}$?" Just use the obtained $a$ as b.c. of the far-field asymptotic equation.

"I tried the following direct shooting method… It gives some errors… " I guess it's because bc1 and bc2 are not b.c in usual meaning, after all, they're equal to 

y[xMax] == (y[xMax] /. DSolve[#, y, x]) & /@ {bc1, bc2}
(* {y[6] == {E^(4/5) C[1]}, y[6] == {C[1] Cos[2/5] + C[2] Sin[2/5]}} *)
Improve this answer
$\endgroup$
13
  • $\begingroup$ @ xzczd, I just noticed your reply. I will learn something from you and reply you later. Many thanks! $\endgroup$
    – W. Robin
    Commented Jan 16, 2016 at 16:49
  • $\begingroup$ @ xzczd, 3 questions, I'm a beginner. 1. Do you remember what value of the $c$ and $xMax$ were used in your answer? 2. I found that you have Map a pure function If[y["Domain"][[1, -1]] < xMax, 1, #] == 0 & to each element in the 1st level of Subtract @@@ {bc1, bc2}. However, applying Subtract to {bc1,bc2} by @@@ is just equivalent to {bc1,bc2}, if I didn't make a mistake on the computation sequence? That is, you apply True or False to {bc1,bc2}? 3. What is the meaning of # in If expression, I know it is commonly used the formal para. for a pure function. Thanks again! $\endgroup$
    – W. Robin
    Commented Jan 17, 2016 at 3:48
  • 1
    $\begingroup$ @W.Robin 1. c = 1; xMax = 6; 2. Observe Subtract @@@ {xx == yy, zz == ww} 3. Yeah, # owns no other meaning, it's the argument of pure function. If you still have trouble in understanding its usage here, just execute If[y["Domain"][[1, -1]] < xMax, 1, #] == 0 & /@ Subtract @@@ {bc1, bc2} outside of ContourPlot and observe the result. $\endgroup$
    – xzczd
    Commented Jan 17, 2016 at 3:58
  • $\begingroup$ @ xzczd, now I understand a little bit: using Subtract @@@ {bc1, bc2} just extract the LHS of each BC. because the 1st level of bc1 or bc2 is just Equal in which the 2nd argument is 0! Am I right? $\endgroup$
    – W. Robin
    Commented Jan 17, 2016 at 4:31
  • 1
    $\begingroup$ @W.Robin The warning isn't a big deal, it's simply because when y["Domain"][[1, -1]] execute for the first time, y["Domain"] isn't yet a list, in v10 a new function Indexed is introduced to handle this problem more properly, but since I (and you) are still in v9, so I simply used Quiet. This doesn't influence the result at all. $\endgroup$
    – xzczd
    Commented Jan 17, 2016 at 5:56

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.