0
$\begingroup$

I have a system of equations to be solved numerically- $$ v'(x) = F_1(v,w,x,\lambda_2,P) \\ w'(x) = F_2(v,w,x,\lambda_2,P) \\ \lambda_2'(x) = \frac{v}{x} $$ The boundary conditions of the problem are-

At $x=0$, $v=w=0$ and $\lambda_2=\lambda_0$ where $\lambda_0$ is some constant.

End the iteration when $\lambda_2=1$ and check value of $x$.

I have found out the system of equations symbolically. I am attaching the code-

    lambda1 = Sqrt[(lambda2[x] + x lambda2'[x])^2 + (eta'[x])^2];
    lambda3 = 1/(lambda1 lambda2[x]);
    I1 = lambda1^2 + lambda2[x]^2 + lambda3^2;
    I2 = 1/lambda1^2 + 1/lambda2[x]^2 + 1/lambda3^2;
    alpha = 0;
    W = -(I1 - 3) - alpha (I2 - 3) + 0.5 x lambda2[x]^2 eta'[x] P;
    eq1 = Numerator[Together[D[W, lambda2[x]] - D[D[W, lambda2'[x]], x]]];
    eq2 = Numerator[Together[D[D[W, eta'[x]], x]]];
    eq1mod = Expand[
       eq1 /. {eta'[x] -> w, eta''[x] -> w', lambda2'[x] -> v/x, 
         lambda2''[x] -> (v' - v/x)/x, lambda2[x] -> lambda2}];
    eq2mod = Numerator[
       Together[
        Expand[eq2 /. {eta'[x] -> w, eta''[x] -> w', lambda2'[x] -> v/x, 
           lambda2''[x] -> (v' - v/x)/x, lambda2[x] -> lambda2}]]];
    f1 = Coefficient[eq1mod, v'];
    f2 = Coefficient[eq1mod, w'];
    f3 = Simplify[-(eq1mod - f1 v' - f2 w')];
    f4 = Coefficient[eq2mod, v'];
    f5 = Coefficient[eq2mod, w'];
    f6 = Simplify[-(eq2mod  - f4 v' - f5 w')];
    F1sym = -(f2 f6 - f3 f5)/(f1 f5 - f2 f4);
    F2sym = (f1 f6 - f3 f4)/(f1 f5 - f2 f4);

I know I am correct upto this point. Now I converted the 2 symbolic expressions to functions F1(x,v,w,lambda2,P) and F2(x,v,w,lambda2,P).

    exprToFunction[expr_, vars_] := 
      ToExpression[
       ToString[FullForm[expr] /. MapIndexed[#1 -> Slot @@ #2 &, vars]] <>
         "&"];
    F1 = exprToFunction[F1sym, {x, v, w, lambda2, P}];
    F2 = exprToFunction[F2sym, {x, v, w, lambda2, P}];

On closer look, I found the system to be singular at x=0. So I started the calculation from x=0.001. I am putting the code below-

    P = 0.5;
    lambda0 = 1.5;
    sol = NDSolve[{
        v'[x] == F1[x, v[x], w[x], lambda2[x], P],
        w'[x] == F2[x, v[x], w[x], lambda2[x], P],
        lambda2'[x] == v[x]/x,
        v[0.001] == 0, w[0.001] == 0, lambda2[0.001] == lambda0,
        WhenEvent[lambda2[x] == 1, "StopIntegration"]},
       {v, w, lambda2}, {x, 0.001, 20}];
    end = InterpolatingFunctionDomain[First[lambda2 /. sol]][[1]]
    Plot[Evaluate[lambda2[x] /. sol], {x, 0.001, 1}, Frame -> True]

Now, after solving, I find the output plot not correct. I don't know the real reason, but I think, in the first step itself, there is a wrong step-size selected. How to correct this problem? The correct image should be this-

Actual Solution

Ignore the multiple plots, because of the parameter, $P$, but the shape should match. I am getting a plot like this-

My_Solution

$\endgroup$
7
  • $\begingroup$ The initial conditions matter a lot. Try setting the start value to where you specify the initial conditions to: 0.001, 0.0009,0.0008, 0.0007, .. and you will see that the system is sensitive to the initial conditions $\endgroup$ Oct 12 '20 at 15:56
  • $\begingroup$ That was a good advice. I tried that, but of no use. Also, I should have specified the initial conditions in the question itself. I am editing the question rn. $\endgroup$
    – user75147
    Oct 12 '20 at 18:13
  • $\begingroup$ I think I know what you mean, but your question is horribly formulated. If you first make a calculation with a start value of 0.1, then 0.01, then 0.001 e.t.c you will see that the solution diverges. Therefore, you can NOT start with 0. $\endgroup$ Oct 12 '20 at 18:44
  • $\begingroup$ Yes, I agree. I could have put up the problem a bit more clearer. As per your advice, I am trying to see the solution by slowly increasing the initial conditions. $\endgroup$
    – user75147
    Oct 12 '20 at 19:47
  • 2
    $\begingroup$ " it should be the step-size, I think." I don't think so. The ODE solver of NDSolve is quite robust and should always be the last thing to tackle. Better to double check the system itself. $\endgroup$
    – xzczd
    Oct 13 '20 at 2:27
2
$\begingroup$

The problem is ill-posed; $v$ cannot vanish at $x = 0$.

Assuming that $v(0) = 0$, $w(0) = 0$, and $\lambda_2(0) = \lambda_0 = 1.5$, we have $$ v'(x) \sim \frac{1.08306}{x} $$ for small $x$, which is inconsistent with $v(0) = 0$.

derivativeList = {
  F1[x, v[x], w[x], lambda2[x], P],
  F2[x, v[x], w[x], lambda2[x], P],
  v[x]/x
};

initialRules = {
  v[x] -> 0,
  w[x] -> 0,
  lambda2[x] -> lambda0
};

derivativeList /. initialRules // FullSimplify
(* {0. + 1.08306/x, 0.308318, 0} *)
$\endgroup$
1
  • $\begingroup$ Thanks for the information. The system of equations coming out has a problem. I have to rethink it. $\endgroup$
    – user75147
    Oct 18 '20 at 4:47

Your Answer

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