3
$\begingroup$

$$\begin{cases}&-z^{\prime\prime}(t)=\lambda(1+(N-2)t)^{\frac1{2-N}(2(N-1)+\alpha)}f(z(t)),\quad t\in(0,+\infty)\\&z(0)=z^\prime(+\infty)=0\end{cases}$$

I'm trying to solve the above differential equation. I'm using the shooting method, but for each initial guess a different solution appears. Which solution is correct? Is it possible to monitor the convergence rate and the iteration steps of the shooting method? Any help or hint is welcome!

Here is the code:

NN = 3.;
p = 3.;
xf = 10000.;

sols = 
  Map[{z -> 
      NDSolveValue[{-z''[
           t] == (1 + (NN - 2) t)^(1/(2 - NN) (2 (NN - 1) + 1)) z[t]^
           p, z[0] == 0, z'[xf] == 0}, z, {t, 0, xf}, 
       Method -> {"Shooting", 
         "StartingInitialConditions" -> {z[0] == 0, z'[0] == #}}]} &,
   {0, 5, 10, 20, 40, 60, 140}];

Plot[Evaluate[z[t] /. sols], {t, 0, 100}, PlotRange -> All]
$\endgroup$
2
$\begingroup$

To assess the correctness of solutions, let's look at how well they meet the boundary conditions you impose.

derivativeAtZero = {0, 5, 10, 20, 40, 60, 140};

With[{NN = 3., p = 3., xf = 10000.}, 
  sols = 
      NDSolveValue[
        {-z''[t] == (1 + (NN - 2) t)^(1/(2 - NN) (2 (NN - 1) + 1)) z[t]^p,
         z[0] == 0, z'[xf] == 0},
        z, {t, 0, xf},
        Method ->
          {"Shooting",
           "StartingInitialConditions" -> {z[0] == 0, z'[0] == #}}] & 
      /@ derivativeAtZero];

computedDerivative = #'[0] & /@ sols;

TableForm[
  Transpose[
    {derivativeAtZero, computedDerivative, derivativeAtZero - computedDerivative}], 
  TableHeadings -> {None, {"Imposed", "Solved", "Difference"}}]

table

The 1st solution, which is the trivial solution z[t] = 0, looks good. For the rest, it's up to you decide if the interpolation function returned by NDSolveValue approximates the imposed boundary condition well-enough to satisfy you. The 2nd, 4th and 5th solutions don't look too bad, so let's plot them.

Plot[Evaluate[#[t] & /@ sols[[{2, 4, 5}]]], {t, 0, 20}, 
  PlotRange -> All, PlotLegends -> {2, 4, 5}]

plot

On general grounds, only the plot corresponding to the 2nd boundary condition looks well behaved, but I can't be sure. I don't know the system being modeled, so I don't really know what kind of transient behavior can be expected.

$\endgroup$
3
$\begingroup$

There is no reason to suppose that the nonlinear ODE presented in the question has a unique solution. Indeed, three were identified here. Several more can be found with

NN = 3; p = 3; xf = 100000; 
sols = Map[{z -> 
    NDSolveValue[{-z''[t] == (1 + (NN - 2) t)^(1/(2 - NN) (2 (NN - 1) + 1)) z[t]^p, 
    z[0] == 0, z'[xf] == 0}, z, {t, 0, xf}, 
    Method -> {"Shooting", "StartingInitialConditions" -> {z[0] == 0, z'[0] == #}}]} &,
   Range[0, 200, 5]];

and then identify the unique solutions.

DeleteDuplicates[sols, Abs[(z'[0] /. #1) - (z'[0] /. #2)] < .01 &];
Select[%, Evaluate[z'[0] /. #] > 0 &];
Sort[%, Evaluate[z'[0] /. #1] < Evaluate[z'[0] /. #2] &];
z'[0] /. %

(* {5.2103, 19.1366, 41.7079, 72.9223, 112.779, 161.279, 218.422, 284.207} *)

LogLinearPlot[Evaluate[z[t] /. %%], {t, 0.01, 1000}, PlotRange -> All]

enter image description here

Solutions have progressively more oscillations before becoming constant. In effect, we have a set of eigenfunctions.

It may be of interest to look at initial-value solutions for a Range of z'[0] near one of these eigenfunctions, for instance,

solt = Map[{z -> 
    NDSolveValue[{-z''[t] == (1 + (NN - 2) t)^(1/(2 - NN) (2 (NN - 1) + 1)) z[t]^p, 
    z[0] == 0, z'[0] == #}, z, {t, 0, 100000}]}&, Range[41.7075, 41.7085, .0001]];
LogLinearPlot[Evaluate[z[t] /. solt], {t, 0.01, xf}, PlotRange -> All]

enter image description here

demonstrating the sensitivity of the asymptotic solution to z'[0].

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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