Skip to main content

differential equations - NDSolve with boundary condition at infinty


I have a nonlinear set of equations with a boundary condition at infinity. Consequently I have to shoot for the boundary. This question has some good ideas but there is an extra complication in my case; my boundary conditon at infinity is the second derivative going to zero and this is approached asymptotically. I will have more problems the same and this is minimum example I can demonstrate.


The equations and initial conditions at zero are here


eqns = {-Derivative[3][a[0]][η] == (-(1/2))*α*
           (-1 + 2*Derivative[1][a[0]][η]^2 + Derivative[1][a[1]][η]^2 + 
              Derivative[1][a[2]][η]^2 + Derivative[1][b[1]][η]^2 + 
              Derivative[1][b[2]][η]^2 - 
       2*a[0][η]*Derivative[2][a[0]][η] - 
              a[1][η]*Derivative[2][a[1]][η] - a[2][η]*Derivative[2][a[2]][

                  η] - b[1][η]*Derivative[2][b[1]][η] - 
              b[2][η]*Derivative[2][b[2]][η]), 
       Derivative[1][b[1]][η] - Derivative[3][a[1]][η] == 
         (1/2)*(-4*α*Derivative[1][a[0]][η]*Derivative[1][a[1]][η] - 
              2*α*Derivative[1][a[1]][η]*Derivative[1][a[2]][η] - 
              2*α*Derivative[1][b[1]][η]*Derivative[1][b[2]][η] + 
              2*α*a[1][η]*Derivative[2][a[0]][η] + 
       2*α*a[0][η]*
                Derivative[2][a[1]][η] + α*a[2][η]*
        Derivative[2][a[1]][η] + 

              α*a[1][η]*Derivative[2][a[2]][η] + α*b[2][η]*
                Derivative[2][b[1]][η] + α*b[1][η]*
        Derivative[2][b[2]][η]), 
       1 - Derivative[1][a[1]][η] - Derivative[3][b[1]][η] == 
         (1/2)*(-4*α*Derivative[1][a[0]][η]*Derivative[1][b[1]][η] + 
              2*α*Derivative[1][a[2]][η]*Derivative[1][b[1]][η] - 
              2*α*Derivative[1][a[1]][η]*Derivative[1][b[2]][η] + 
              2*α*b[1][η]*Derivative[2][a[0]][η] + 
              α*b[2][η]*Derivative[2][a[1]][η] - α*b[1][η]*
                Derivative[2][a[2]][η] + 

       2*α*a[0][η]*Derivative[2][b[1]][η] - 
              α*a[2][η]*Derivative[2][b[1]][η] + α*a[1][η]*
                Derivative[2][b[2]][η]), 2*Derivative[1][b[2]][η] - 
           Derivative[3][a[2]][η] == 
         (1/2)*(α - α*Derivative[1][a[1]][η]^2 - 
              4*α*Derivative[1][a[0]][η]*Derivative[1][a[2]][η] + 
              α*Derivative[1][b[1]][η]^2 + 
       2*α*a[2][η]*Derivative[2][a[0]][
                  η] + α*a[1][η]*Derivative[2][a[1]][η] + 
              2*α*a[0][η]*Derivative[2][a[2]][η] - 

              α*b[1][η]*Derivative[2][b[1]][η]), 
       -2*Derivative[1][a[2]][η] - Derivative[3][b[2]][η] == 
         (1/2)*(-2*α*Derivative[1][a[1]][η]*Derivative[1][b[1]][η] - 
              4*α*Derivative[1][a[0]][η]*Derivative[1][b[2]][η] + 
              2*α*b[2][η]*Derivative[2][a[0]][η] + 
              α*b[1][η]*Derivative[2][a[1]][η] + α*a[1][η]*
                Derivative[2][b[1]][η] + 2*α*a[0][η]*Derivative[2][b[2]][
                  η])}; 
ic0 = {a[0][0] == 0, Derivative[1][a[0]][0] == 0, a[1][0] == 0, 
       b[1][0] == 0, Derivative[1][a[1]][0] == 0, Derivative[1][b[1]][0] == 

         0, a[2][0] == 0, b[2][0] == 0, Derivative[1][a[2]][0] == 0, 
       Derivative[1][b[2]][0] == 0};

I then have a shooting function that determines the values of the second derivative at a point that is getting towards infinity.


ClearAll[shoot]; 
shoot[{(a0_)?NumberQ, a1_, b1_, a2_, b2_}, inf_] := 
   (sol = First[NDSolve[Evaluate[Join[eqns, ic0, 
                 {DerivatClearAll[shoot]; 
shoot[{(a0_)?NumberQ, a1_, b1_, a2_, b2_}, inf_] := 
   (sol = First[NDSolve[Evaluate[Join[eqns, ic0, 

                 {Derivative[2][a[0]][0] == a0, Derivative[2][a[1]][0] == a1, 
                   Derivative[2][b[1]][0] == b1, Derivative[2][a[2]][0] == a2, 
                   Derivative[2][b[2]][0] == b2}] /. {α -> 0.1}], 
           {a[0], a[1], b[1], a[2], b[2], Derivative[2][a[0]], 
             Derivative[2][a[1]], Derivative[2][b[1]], Derivative[2][a[2]], 
             Derivative[2][b[2]]}, {η, 0, inf}]]; 
     res = {Derivative[2][a[0]], Derivative[2][a[1]], Derivative[2][b[1]], 
           Derivative[2][a[2]], Derivative[2][b[2]]} /. sol; 
     (#1[inf] & ) /@ res)


By playing around I have found some good starting values and when I use these with infinity set to 9 everything works.


starts = {{a0, 0.054}, {a1, 0.704}, {b1, -0.702}, {a2, 0.021}, {b2, 
    0.020}};
inf = 9;
rts = FindRoot[shoot[{a0, a1, b1, a2, b2}, inf] == {0, 0, 0, 0, 0}, 
   starts];
shoot[{a0, a1, b1, a2, b2} /. rts, inf]


{-6.60792*10^-14, 5.00619*10^-13, 3.11087*10^-13, -3.76011*10^-13, 1.65523*10^-13}




If I plot the second derivative I can see that it goes to zero.


Join[{Plot[Evaluate[{a[0]''[η]} /. sol], {η, 0, inf}, 
   PlotRange -> All]},
 Table[
  Plot[Evaluate[{a[n]''[η], b[n]''[η]} /. sol], {η, 0, 
    inf}, PlotRange -> All],
  {n, 1, 2}
  ]
 ]


Mathematica graphics


However my infinity is not big enough. If I increase it to say 10 or more I get


inf = 10;
rts = FindRoot[shoot[{a0, a1, b1, a2, b2}, inf] == {0, 0, 0, 0, 0}, 
   starts];
shoot[{a0, a1, b1, a2, b2} /. rts, inf]


FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.




If one runs the integration on from the results with inf = 9 then one can see that the second derivative was made to pass through the point 9 but is not zero subsequently. Here I have run the inf = 9 case to get the roots but now solve and plot using these values for inf = 15.


inf = 15;
sol = First@NDSolve[Evaluate[Join[eqns, ic0,
       {(a[0]'')[0] == 
         a0, (a[1]'')[0] == 
         a1, (b[1]'')[0] == b1,
        (a[2]'')[0] == 
         a2, (b[2]'')[0] == b2}] /. 
      Join[{α -> 0.1}, rts]],

    {a[0], a[1], b[1], a[2], b[2], a[0]'', a[1]'', b[1]'', a[2]'', 
     b[2]''}, {η, 0, inf}];
Join[{Plot[Evaluate[{a[0]''[η]} /. sol], {η, 0, inf}, 
   PlotRange -> All]},
 Table[
  Plot[Evaluate[{a[n]''[η], b[n]''[η]} /. sol], {η, 0, 
    inf}, PlotRange -> All],
  {n, 1, 2}
  ]
 ]


Mathematica graphics


So I have failed to get the asymptotic solution. I have tried integrating over an interval at large values of η and then using FindMinimum but this does not give good results.


Do you have any further ideas? Thanks



Answer



Once again, let's turn to the primitive finite difference method (FDM). I'll use pdetoae for the generation of difference equation:


var = Flatten@{a /@ Range[0, 2], b /@ Range@2}
bc = D[Through@var[t] == 0 // Thread, t, t] /. t -> inf

inf = 15;

points = 200;
difforder = 4;
grid = Array[# &, points, {0, inf}];
ptoafunc = pdetoae[var[η], grid, difforder];
del = #[[3 ;; -2]] &;
ae = del /@ ptoafunc@eqns;
(* Definition of pdetoae isn't included in this post,
   please find it in the link above. *)
aebc = ptoafunc@{ic0, bc};


sollst = Partition[
    With[{guess = 1}, 
      FindRoot[{ae, aebc} /. {α -> 1/10}, 
       Flatten[Table[{var[[i]][η], guess}, {i, 5}, {η, grid}], 1]]][[All, -1]],
     points]; // AbsoluteTiming


Remark


The following way to calculate sollst doesn't require one to use del:


fullsys = Flatten@ptoafunc@{eqns, ic0, bc};


lSSolve[obj_List, constr___, x_, opt : OptionsPattern[FindMinimum]] := 
 FindMinimum[{1/2 obj^2 // Total, constr}, x, opt]
lSSolve[obj_, rest__] := lSSolve[{obj}, rest]

sollst = Partition[
    With[{guess = 1}, 
      lSSolve[Subtract @@@ fullsys /. {α -> 1/10}, 
       Flatten[Table[{var[[i]][η], guess}, {i, 5}, {η, grid}], 1]]][[2, 
      All, -1]], points]; // AbsoluteTiming


But notice it's slower and less robut compared to the FindRoot approach.



solfunclst = ListInterpolation[#, grid] & /@ sollst;

Plot[D[Through@solfunclst[x], x, x] // Evaluate, {x, 0, inf}, PlotRange -> All]

Mathematica graphics


FindRoot still spits out lstol, but further check by adjusting points, difforder, etc. shows the result seems to be accurate enough.


Comments

Post a Comment

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...
Post a Comment
Read more

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...
Post a Comment
Read more