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.

With this procedure, one may determine an eigen-value function $R(a)$ for any given $\Xi$ (say 0, 25, 50, 75, 100)

g[eta1_, eta3_, R_, a_, Xi_] := {f[1], h[1], q[1]} /.     
NDSolve[{D[f[x], {x, 2}] - a^2 f[x] + a (h[x] + R q[x]) == 0,
    D[h[x], { x, 2}] - a^2 h[x] + a Xi x f[x] == 0,
    D[q[x], {x, 2}] - a^2 q[x] + a f[x] == 0,
     f[0] == 0,
f'[0] == 1,
h[0] == eta1,
h'[0] == 0, 
q[0] == 0,
q'[0] == eta3 
}, {f, h, q}, {x, 0, 1}]
s[a_, Xi_] :=FindRoot[g[eta1, eta3, R, a, Xi], {{R, 0}, {eta1, 0}, {eta3, 0}}, 
  MaxIterations -> 100000]

where the condition $f'(0) = 1$ is a normalisation condition. But the problem is that the code doesn't converge always. The output I am looking for is like this

enter image description here

Thanks

share|improve this question
1  
Please add some code to demonstrate where your implementation allegedly fails. –  Yves Klett Oct 28 '13 at 13:57
    
Well, for s[7.8,100], R=28 and for s[8,100], R=22, which is not what we can see from the above fig. –  MMM Oct 28 '13 at 14:26
    
Please update your question accordingly and give us a simple template to check against the figure. –  Yves Klett Oct 28 '13 at 14:46
    
I noticed my mistake and rewrited my answer, I think now it's just what you want. Have a look! –  xzczd Nov 4 '13 at 13:11
add comment

1 Answer 1

up vote 4 down vote accepted
+50

Except for simple cases (the only case I can think of is monotonic function at the moment), results of FindRoot are often affected by the starting point for root searching. More starting points, more uncertain, so FindRoot is easier to fail when it's searching for the roots of a set of equations than that of a single equation. The purpose of your code is in fact searching $R$ that makes $f(1)$, $h(1)$, $q(1)$ equal to $0$ for given $a$ and $Ξ$, i.e. $f(x)$, $h(x)$, $q(x)$ will satisfy the boundary condition(B.C.):

$$f(1)=h(1)=q(1)=0$$

Why not use some of these B.C.s instead of the uncertain ones inside NDSolve?

After some trial I found that using $f(1)=q(1)=0$ inside seems to be the best choice, so I modified your code into:

(*
 h[0] == eta1 => f[1] == 0
q'[0] == eta3 => q[1] == 0
*)
g2[R_?NumericQ, a_?NumericQ, Xi_?NumericQ] := 
 h[1] /. NDSolve[{D[f[x], {x, 2}] - a^2 f[x] + a (h[x] + R q[x]) == 0,
                  D[h[x], {x, 2}] - a^2 h[x] + a Xi x f[x] == 0, 
                  D[q[x], {x, 2}] - a^2 q[x] + a f[x] == 0, 
                  f'[0] == 1, h'[0] == 0, q[0] == 0, f[0] == 0, f[1] == 0, q[1] == 0}, 
                 {f, h, q}, {x, 0, 1}]
s2[a_, Xi_] := FindRoot[g2[R, a, Xi], {R, 0}]

The role of those ?NumericQs is to quit the warning NDSolve::ndinnt and ReplaceAll::reps, the output won't be affected even if you take them away.

Let's check the effect of the new definition:

sample = Join[Range[6/10, 2, 1/10], Range[2, 8, 1/4]];
data[Xi_] := {sample, R /. s2[#, Xi] & /@ sample} // Transpose;
rst = data /@ Range[0, 100, 25];
ListLinePlot[rst, PlotRange -> {{0, 8}, {-40, 100}}, Frame -> True, Axes -> False]

Some warnings still generate, but the result matches your figure very well:

enter image description here

And needless to say, much better than your original:

olddata[Xi_] := {sample, R /. s[#, Xi] & /@ sample} // Transpose;
oldrst = olddata /@ Range[0, 100, 25];
ListLinePlot[oldrst, PlotRange -> {{0, 8}, {-40, 100}}, Frame -> True, Axes -> False]

enter image description here

share|improve this answer
    
lovely! thx mate –  MMM Nov 4 '13 at 17:05
add comment

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.