4

I'm trying to solve a pair of coupled ODE's. I need to place four Dirichlet boundary conditions (at R = 0 and R = ∞ for each dependent variable), and my independent variable R runs over R > 0. I'm approximating these boundary conditions with a large and a small value of R to avoid a singularity at R = 0. I know the system has a well behaved, smooth solution for X and P. (In fact X is monotonically increasing and P is monotonically decreasing and both should be virtually constant by about R = 10.)

When I run NDSolve however, Mathematica throws

NDSolve::ndsz: At R == 0.034300651057220126`, step size is effectively zero; singularity or stiff system suspected.

I don't believe this system should be stiff, and there shouldn't be any singularities in the range I'm evaluating over.

I have tried modifying the working precision, the number of steps, the method, the starting step size, all to no avail.

Any help would be appreciated.

eqn1 = -X''[R] - X'[R]/R + (P[R]^2 X[R])/R^2 + 1/2 X[R] (X[R]^2 - 1) == 0
eqn2 = -P''[R] + P'[R]/R + (X[R]^2 P[R]) == 0


inf = 100; 
zero = 0.0001; 
eqns := {eqn1, eqn2}; 
bcs = { X[inf] == 1,
        P[inf] == 0,
        X[zero] == 0,
        P[zero] == 1}


sol = NDSolve[eqns~Join~bcs, {X, P}, {R, zero, inf}, Method -> Automatic]


p1 = Plot[Evaluate[X[R] /. sol], {R, zero, 10}, PlotRange -> All];
p2 = Plot[Evaluate[P[R] /. sol], {R, zero, 10}, PlotRange -> All];
Show[p1, p2]
A. Wills.
  • 43
  • 4
  • I do love a good topological defect problem. The problem stems from the fact that Mathematica uses shooting methods for BVPs, but the solutions for equations of this type for generic initial conditions often do diverge. Let me see if I can dig up some guidance on this type of problem. – Michael Seifert May 25 '18 at 17:44
  • Ah, here's the answer I was thinking of: Numerical solution of coupled ODEs with boundary conditions, particularly bbgodfrey's answer there. Note the lengths to which he had to go to cajole Mathematica into solving this problem. That said, on a practical level, my advice to you is to use MATLAB if you're familiar with it; it has several relaxation solvers that allow for non-linear ODEs. (Mathematica's only allow for linear systems.) – Michael Seifert May 25 '18 at 18:20
  • @KraZug: Other than this problem not depending on two free parameters, the problems look quite similar to me. But of course they're not identical, and there might be some properties of this system that make it easier to solve. If I'm wrong and you or someone else posts a simpler solution, I'll happily eat crow. – Michael Seifert May 25 '18 at 18:30
  • @MichaelSeifert, I think it is me that needs to eat crow. My apologies, I thought the equations were simpler than they are. – SPPearce May 25 '18 at 18:47
  • As suggested in the comments above, it may take heroic effort to solve these two second-order, nonlinear, non-autonomous ODEs. To get started, derive symbolically the first few terms of a power series solution at r == 0 and the asymptotic solution at large r. The problem then becomes connecting the two. – bbgodfrey May 25 '18 at 22:55
  • In fact, the equations here are very similar to those in my solution mentioned by @MichaelSeifert above, the only significant difference being that ξ^2 must be replaced by -1. Boundary conditions also are different, but probably not significantly so. – bbgodfrey May 26 '18 at 01:07
  • I note that the second terms in each equation have differing signs, and that the third term in the first equation is divided by $R^2$ while the corresponding term in the second equation is not. Is this correct? – Michael Seifert May 26 '18 at 14:23
  • Michael - Thank you for your help! The system in the post you linked is indeed the same as the one I am attempting to solve up to variable redefinitions, with slightly different boundary conditions. There was a sign error in the third term of eqn2 which I have now edited in the post. You now will probably recognise the problem, since you spotted its origin! Indeed the solution is known, this was more an exercise for myself in using Mathematica than anything else.

    bbgodfrey - Thank you for your suggestions and the thorough answer in the linked question. I believe it will help immensely!

     – A. Wills. May 26 '18 at 15:01

1 Answers1

2

After the equation is corrected, FDM approach works, so let's step backward once again. I'll use pdetoae for the generation of difference equation:

eqn1 = -X''[R] - X'[R]/R + (P[R]^2 X[R])/R^2 + 1/2 X[R] (X[R]^2 - 1) == 0;
eqn2 = -P''[R] + P'[R]/R + (X[R]^2 P[R]) == 0;
inf = 20;
zero = 0;
eqns = Map[R^2 # &, {eqn1, eqn2}, {2}] // Simplify;
bcs = {X[inf] == 1, P[inf] == 0, X[zero] == 0, P[zero] == 1};
domain = {zero, inf};
points = 100;
grid = Array[# &, points, domain];
difforder = 2;
(* Definition of pdetoae isn't included in this post,
   please find it in the link above. *)
ptoafunc = pdetoae[{X, P}[R], grid, difforder];
del = #[[2 ;; -2]] &;

ae = del /@ ptoafunc@eqns;

initialguess = 1;        
solrule = FindRoot[{ae, bcs}, 
   Flatten[#, 1] &@Table[{var[R], initialguess}, {var, {X, P}}, {R, grid}]];

{sol@X, sol@P} = ListInterpolation[#, domain] & /@ Partition[solrule[[All, -1]], points];

Plot[{sol[X][R], sol[P][R]}, {R, zero, inf}]

Mathematica graphics

xzczd
  • 65,995
  • 9
  • 163
  • 468
  • This program does not give any GRAPH. Actually I am new to MATHEMATICA, After copy and paste this code and running the code, no graph is coming I dont know why? Please help – smohanty May 02 '19 at 03:46
  • @smohanty As already mentioned above, definition of pdetoae isn't included in this answer, please find it in the link. – xzczd May 02 '19 at 05:18