Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

my equation is:

k[a_?NumericQ, b_?NumericQ, c_?NumericQ] := NDSolve[{r^3 b x''[r] y[r] - 12 x[r] y[r] - 
 4 r^3 b y''[r] x[r] == 0, 4 (2 x'[r]^2 + 
    2 x[r]x''[r])/Sqrt[c] + (2 y'[r]^2 + 
    2 y[r] y''[r])/Sqrt[c] == -4b/r,x[1] == 0, y[1] == Sqrt[-4 b Log[2/b]/(c^(-1/2))],
y'[1] == c^(1/4) (-4 b - 4 b Log[2 /b])/(2 Sqrt[-4 b Log[2/b] ] ), x'[1] == a},
 {x, y}, {r, 1, b/2}, MaxSteps -> 1000000]

and then I plot it:

Plot[Evaluate[{x[r], y[r]} /. k[1, 50, 10^(-7)]], {r, 1, 25}]

enter image description here

this one is okay, cause x[r] and y[r] both equal to 0 at r=b/2

if I set a=0 and the other parameters stay the same

enter image description here

I do not know how comes this sigularity at r=b/2, basing on the second equation in NDSolve, 4 (2 x'[r]^2 + 2 x[r]x''[r])/Sqrt[c] + (2 y'[r]^2 + 2 y[r] y''[r])/Sqrt[c] == -4b/r comes from the second derivative of r in this equation:

c^(-1/2) 4 x[r]^2 + c^(-1/2) y[r]^2 == -4 r b Log[(2 r)/b]

so it's easy to check, if x[r]=0 all the way, y[r] should be 0 at r=b/2

does anyone know how to resolve the infinity at r=b/2 for any values of the parameter? a>=0, b>2, 10^(-3)>c>10^(-7)

thanks a lot

share|improve this question

1 Answer 1

up vote 2 down vote accepted

I think this is just another story of WorkingPrecision, let's try a higher one, say, 32:

k[a_?NumericQ, b_?NumericQ, c_?NumericQ] := 
 NDSolve[{r^3 b x''[r] y[r] - 12 x[r] y[r] - 4 r^3 b y''[r] x[r] == 0,
    4 (2 x'[r]^2 + 2 x[r] x''[r])/
       Sqrt[c] + (2 y'[r]^2 + 2 y[r] y''[r])/Sqrt[c] == -4 b/r, 
   x[1] == 0, y[1] == Sqrt[-4 b Log[2/b]/(c^(-1/2))], 
   y'[1] == c^(1/4) (-4 b - 4 b Log[2/b])/(2 Sqrt[-4 b Log[2/b]]), 
   x'[1] == a}, {x, y}, {r, 1, b/2}, WorkingPrecision -> 32]

Plot[Evaluate[{x[r], y[r]} /. k[0, 50, 10^(-7)]], {r, 1, 25}]

enter image description here

OK, the singularity disappears.

share|improve this answer
    
yes, it works. can you explain more to me why this works? thanks –  3c. Jul 26 '13 at 14:07
    
@3c. Well, to be honest, I'm not able to tell the exact reason, but based on my experience, if the plot of a function becomes unusual(singularity, hard oscillation etc while the function should not possess these features) in a point or in a small interval, it's often related to a WorkingPrecision that isn't high enough. There're numbers of posts about WorkingPrecision in this site, e.g. this, you can have a search for more details. –  xzczd Jul 26 '13 at 14:26
    
thanks very much –  3c. Jul 26 '13 at 14:36

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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