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.

Here is my code:

k = 0.42; e = 5; d = 7; n = 1/0.7158; R2 = 40000; 

s = 
 NDSolve[{k^2  (u'[y] )^(3/2) (u'[y] + 2 u[y] Exp[y]/(R2 - Exp[y]))^(
      1/2) - (-u'[y] (R2 Exp[-y] - 1) - 2 u[y]) (n - 1)/R2 - 1 == 0, 
   u[Log[4.17 (n - 1)]] == 4.17}, 
  u, {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}];

C1[y_] = Evaluate[u[y] /. s[[2]]];

Plot[C1[y], {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]

x = (((-C1'[y] (R2 Exp[-y] - 1) - 2 C1[y])/(-Exp[y] + R2)^3 + 
      R2/((n - 1) (-Exp[y] + R2)^3)) Abs[
     C1'[y] Exp[-y] + 2 C1[y]/(R2 - Exp[y])])^(-1/4) Exp[-y] (n - 1)^(
  1/2) R2^(-1/2)

f[y_] = NIntegrate[
  2/3 a^(-5/3)
    Exp[-a d Abs[x]] (1 + (e n Exp[y]/(R2 (n - 1)))^2/a^2)^(-17/
    6), {a, 1, Infinity}]

f1 = FunctionInterpolation[(f[y])^(
  3/4), {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]

Plot[Re[f1[y]], {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]

s = NDSolve[{k^2  (f1[y]) (u'[y])^(
      3/2) (u'[y] + 2 u[y] Exp[y]/(R2 - Exp[y]))^(
      1/2) - (-u'[y] (R2 Exp[-y] - 1) - 2 u[y]) (n - 1)/R2 - 1 == 0, 
   u[Log[4.17 (n - 1)]] == 4.17}, 
  u, {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]

C1[y_] = Evaluate[u[y] /. s[[2]]]

Plot[Re[C1[y]], {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]

x = (((-C1'[y] (R2 Exp[-y] - 1) - 2 C1[y])/(-Exp[y] + R2)^3 + 
      R2/((n - 1) (-Exp[y] + R2)^3)) Abs[
     C1'[y] Exp[-y] + 2 C1[y]/(R2 - Exp[y])])^(-1/4) Exp[-y] (n - 1)^(
  1/2) R2^(-1/2)

f[y_] = NIntegrate[
  2/3 a^(-5/3)
    Exp[-a d Abs[x]] (1 + (e n Exp[y]/(R2 (n - 1)))^2/a^2)^(-17/
    6), {a, 1, Infinity}]

f1 = FunctionInterpolation[(f[y])^(
  3/4), {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]


Plot[Re[f1[y]], {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]

(*up to here, I think there is no problem. then the kernel crashes below*)

s = NDSolve[{k^2  (f1[y]) (u'[y])^(
      3/2) (u'[y] + 2 u[y] Exp[y]/(R2 - Exp[y]))^(
      1/2) - (-u'[y] (R2 Exp[-y] - 1) - 2 u[y]) (n - 1)/R2 - 1 == 0, 
   u[Log[4.17 (n - 1)]] == 4.17}, 
  u, {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]

C1[y_] = Evaluate[u[y] /. s[[2]]]

Plot[Re[C1[y]], {y, Log[4.17 (n - 1)], Log[(n - 1) R2/(2 n)]}]

R2 is a constant. When I set R2 smaller than 33000, I can always got a numerical solution of u[y]. However when R2 gets very large,i.e. 40000, mathematica just doesn't work at all. I mean there is no solution given for the third differential equation.

Can anybody tell me how to solve this problem? I will be more than appreciated.

Thanks.

share|improve this question
    
Please post a fully self contained non-working example that one can copy/paste and run to show the problem including any numerical values used. –  Nasser Mar 11 at 0:57
    
Thanks for your reply. I have posted the full code. –  zeo Mar 11 at 1:07
    
Can you reduce the size of your example, and indicate where exactly and what goes wrong? It'll save people a lot of time and might get you an answer. All I see now is that at some point the kernel crashes, but I didn't take the time to find out at which step. I think you can do this yourself and let us know. –  Szabolcs Mar 11 at 1:32
    
That said, a kernel crash indicates a bug, so can you report this to support at wolfram.com ? –  Szabolcs Mar 11 at 1:34
    
Thanks, I edited the post again, indicating where exactly goes wrong –  zeo Mar 11 at 1:45

Your Answer

 
discard

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

Browse other questions tagged or ask your own question.