0

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
Community
  • 1
ApplMath
  • 71
  • 4
  • The first code you posted works OK here. Try restarting Mathematica. If it doesn't work, post the version of Mathematica you're using – Dr. 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 Answers1

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.

ApplMath
  • 71
  • 4