4
$\begingroup$

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]
|improve this question|||||
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Michael Seifert May 25 '18 at 17:44
  • $\begingroup$ 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.) $\endgroup$ – Michael Seifert May 25 '18 at 18:20
  • $\begingroup$ @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. $\endgroup$ – Michael Seifert May 25 '18 at 18:30
  • $\begingroup$ @MichaelSeifert, I think it is me that needs to eat crow. My apologies, I thought the equations were simpler than they are. $\endgroup$ – KraZug May 25 '18 at 18:47
  • $\begingroup$ 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. $\endgroup$ – bbgodfrey May 25 '18 at 22:55
1
$\begingroup$

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

|improve this answer|||||
$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – smohanty May 2 '19 at 3:46
  • $\begingroup$ @smohanty As already mentioned above, definition of pdetoae isn't included in this answer, please find it in the link. $\endgroup$ – xzczd May 2 '19 at 5:18

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.