9
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ Please add some code to demonstrate where your implementation allegedly fails. $\endgroup$ – Yves Klett Oct 28 '13 at 13:57
  • $\begingroup$ 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. $\endgroup$ – zhk Oct 28 '13 at 14:26
  • $\begingroup$ Please update your question accordingly and give us a simple template to check against the figure. $\endgroup$ – Yves Klett Oct 28 '13 at 14:46
  • $\begingroup$ I noticed my mistake and rewrited my answer, I think now it's just what you want. Have a look! $\endgroup$ – xzczd Nov 4 '13 at 13:11
4
+50
$\begingroup$

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

| improve this answer | |
$\endgroup$
  • $\begingroup$ lovely! thx mate $\endgroup$ – zhk Nov 4 '13 at 17:05

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.