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.

I am trying to use a solution given by Michael E2 in this topic:

Zeros[z_] = (BesselJ[0, z] - BesselJ[2, z])*(3*BesselI[1, z] + 
  BesselI[3, z]) - (BesselI[0, z] + 
  BesselI[2, z])*(3*BesselJ[1, z] - BesselJ[3, z]);

(zeros = First @ Last @ Reap @ Quiet @ NDSolve[
   {y'[x] == Zeros'[x], y[0] == Zeros[0], 
    WhenEvent[y[x] == 0, Sow[FindRoot[Zeros[z], {z, x}]]]}, 
   y, {x, 0, 40}, AccuracyGoal -> 1, PrecisionGoal -> 1]) // AbsoluteTiming

-- to find roots for my equation:

$$ J_1(a/z)\, Y_0(b/z)-J_0(a/z)\,Y_1(b/z)=0, \quad a=10, \;b=1, $$

but get an error

First::first: "{} has a length of zero and no first element."

The problem is that "y(x)=0" event doesn't happen. I understand that the functions are quite different, so I perform the following:

$$ J_1(a/z)\; Y_0(b/z)-J_0(a/z)\; Y_1(b/z)=0 \rightarrow\\ \rightarrow \left( J_1(a\tilde z)\; Y_0(b\tilde z)-J_0(a\tilde z)\; Y_1(b\tilde z) \right)\; e^{-100/\tilde z}=0, ~\tilde z = 1/z $$

to make my function to be more like the original one, but still get the error. Could anyone help, please?

My code (I mean, the original code with my function) is below.

Sigma = 1;
xi1 = 10;
Zeros[z_] = 
  (BesselJ[1, xi1*z]*BesselY[0, Sigma*z] - 
   BesselJ[0, Sigma*z]*BesselY[1, xi1*z])*Exp[-100/z];
(zeros = First@
Last@Reap@
  Quiet@NDSolve[{y'[x] == Zeros'[x], y[0] == Zeros[0], 
 WhenEvent[y[x] == 0, Sow[FindRoot[Zeros[z], {z, x}]]]}, 
y, {x, 0, 10}, AccuracyGoal -> 10, PrecisionGoal -> 10]) // AbsoluteTiming
share|improve this question
    
The first code you posted works OK here. Try restarting Mathematica. If it doesn't work, post the version of Mathematica you're using –  belisarius Nov 20 '13 at 20:03
    
@belisarius First code works good. The problem is that I can't get it working with another function, even if I make that function look like the first one (more or less). See second code. My version of Mathematica is 9.0.1.0, 64-bit. –  ApplMath Nov 23 '13 at 10:51

1 Answer 1

Solved it! Here's the code.

xi1 = 10;
Sigma = 1;
f3[z_] = BesselJ[1, xi1/z]*BesselY[0, Sigma/z] - 
   BesselJ[0, Sigma/z]*BesselY[1, xi1/z];

Timing[roots3 = 
Reap[First[
 NDSolve[{x'[z] == f3'[z], x[1] == f3[1]}, {x}, {z, 15, 0.1}, 
  Method -> {"EventLocator", "Event" -> x[z], 
    "EventAction" :> Sow[z]}, MaxStepSize -> 0.001, 
  MaxSteps -> 1000000]]][[2, 1]]]

And

f3 /@ roots3

shows that the roots are quite accurate. Sure it doesn't find all the roots in the interval and skips some of them, but one can easily increase the accuracy.

Note: I used my initial function.

share|improve this answer

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.