NDSolve inside For loop, reseting conditions after each zero

Questions
1

Hi guys, basically what I am trying to do is see whether my graphical solution to an ODE will give a zero at the correct point. For a certain set of initial conditions, this should be the case up until some point where the zeros will vanish (or occur after that specified point). I tried to do this on code but it is not working how I thought - I want to have it see whether x[t-spec] == 0 and if it does, then I will add a certain value (0.000015) to “r” (from the code) and check the next set of initial conditions (and if they satisfy the x[t-spec] == 0 as well, then add a 0.000015 to r). I am interested in finding the “r” value at which beyond it, there will be no more zeros at that specified point. Here is my code below:

    r = 0 (*defines the critical radius variable*)
For[i = 0.000015, i = 1.5, i += 0.000015] { (*iterates through initial condition values*)
  s = NDSolve[{x'[
         t] - (2.5*10^-25) ((1.2*10^31 ((9*10^-6) - (1.5*10^-3) t) \
(5*10^5 ((9*10^-6) - (1.5*10^-3) t) - 1000000 x[t] + 
               3) (5*10^5 ((9*10^-6) - (1.5*10^-3) t) + 
                1000000 x[t] - 
                3)/(1000000000000 ((9*10^-6) - (1.5*10^-3) t)^2 + 
                 1000000000000 x[t]^2 - 6000000 x[t] + 9)^(7/
                 2)) - ((3 ((9*10^-6) - (1.5*10^-3) t)^3 - 
               12 ((9*10^-6) - (1.5*10^-3) t) x[t]^2))/(x[
                 t]^2 + ((9*10^-6) - (1.5*10^-3) t)^2)^(7/2)) == 0, 
      x[0] == i*10^-6}, x, {t, 0, 1}]
    Plot[Evaluate[x[t] /. s], {t, 0, .1}]    (*plots sol'n to diff eq*)
    If[{x[.006] /. s} == 0, r = r + 0.000015, r = r] (*checks if zero at specified time - adds .000015 to r if true)
  }
Print[r]

The logic makes sense but it seems like I am having trouble with the syntax. I am most concerned with the last two lines in the loop - where I plot and evaluate to see whether there is a zero at the specified time. Any help with this would be needed! Also, whenever I run it, I keep getting it showing up as r = 0, so I’m guessing it’s not adding to “r” - that’s the only clue I can see so far. Thanks for the help! Please reply if you have any questions about this so I can clear it up if you’re unsure on something I said!

2

From your example code above, I’m guessing you have had some prior experience in another programming language and are trying to use some of that to write Mathematica.

First problem

For[n = 1, n = 3, n++]

does almost nothing and is done.

This does a little more

For[n = 1, n <= 3, n++]

but not much more.

This does a little more

For[n = 1, n <= 3, n++,
 Print[n];
 Print[n^2]
 ]

Note that Mathematica does not use { and } to group multiple statements, it uses semicolons to accomplish that. I am guessing your prior programming experience with other “more normal” programming languages is leading you to expect to group statements with { and }.

r = 0;(*defines the critical radius variable*)
For[i = 0.000015, i <= 0.00015, i += 0.000015,(*iterates through initial condition values*)
 s = NDSolve[{x'[t] - (2.5*10^-25) ((1.2*10^31 ((9*10^-6) - (1.5*10^-3) t) (5*10^5 ((9*10^-6) - (1.5*10^-3) t) -
    1000000 x[t] + 3)(5*10^5 ((9*10^-6) - (1.5*10^-3) t) + 1000000 x[t] - 3)/(1000000000000 ((9*10^-6) -
    (1.5*10^-3) t)^2 + 1000000000000 x[t]^2 -6000000 x[t] + 9)^(7/2)) - ((3 ((9*10^-6) - (1.5*10^-3) t)^3 -
    12 ((9*10^-6) - (1.5*10^-3) t) x[t]^2))/(x[t]^2 + ((9*10^-6) - (1.5*10^-3) t)^2)^(7/2)) == 0, 
     x[0] == i*10^-6}, x, {t, 0, 1}][[1]];
 Print[Plot[Evaluate[x[t] /. s], {t, 0, .1}]];
 Print["x[.006]=", x[t] /. s /. t -> .006];
 If[x[t] /. s /. t -> .006 == 0, r = r + 0.000015](*checks if zero at specified time-adds .000015 to r if true*)
 ];
Print[r]

I’m not sure you really mean to do 100000 plots so I reduced 1.5 to 0.00015 for a first try.
There are several other changes in that and you should probably compare that with your original, character by character, word by word, and try to figure out why each change was made. But at least it appears to be correctly iterating and displaying plots for you. If I have misunderstood then please let me know and I will try to fix what I can. Please test this very carefully to make certain that the result is correct before depending on it.

I am concerned that your If statement may never find x[t] to be exactly 0 and thus never True. That may not be good programming.

3

Hey Bill, I am having some trouble still. Sorry for the delay.

r = 0;(*defines the critical radius variable*)
For[i = 1.49999, i <= 1.5, 
   i += 0.000015](*iterates through initial condition values*)s = 
 NDSolve[{x'[
      t] - (2.5*10^-25) ((1.2*10^31 ((9*10^-6) - (1.5*10^-3) t) \
(5*10^5 ((9*10^-6) - (1.5*10^-3) t) - 1000000 x[t] + 
            3) (5*10^5 ((9*10^-6) - (1.5*10^-3) t) + 1000000 x[t] - 
             3)/(1000000000000 ((9*10^-6) - (1.5*10^-3) t)^2 + 
              1000000000000 x[t]^2 - 6000000 x[t] + 9)^(7/
              2)) - ((3 ((9*10^-6) - (1.5*10^-3) t)^3 - 
            12 ((9*10^-6) - (1.5*10^-3) t) x[t]^2))/(x[
              t]^2 + ((9*10^-6) - (1.5*10^-3) t)^2)^(7/2)) == 0, 
   x[0] == i*10^-6}, 
  x, {t, 0, 
   1}]; (*defines the differential equation to be solved with init. \
conditions*)
Print["x[.006]=", x[t] /. s /. t -> .006];
If[x[t] /. s /. t -> .006 <= 10^-10, r = r + 0.000015, 
 Print[Plot[
   Evaluate[x[t] /. s], {t, 
    0, .1}]]];(*checks if zero at specified time-adds .000015 to r if \
true*)
Print[r]

The goal of this code is to be able to solve this differential equations for the given set of initial conditions, and produce an ‘r’ value known as the critical radius. Basically, a particle is moving according to this differential equation, and if it reaches within a certain distance (in this case, 10-10 m) of the trap, it will be captured. The differential equation dictates the position of this particle. The range of positions are from 0 to 1.5 units, and what I want to do is find out what the critical radius is – that is, what is the distance ‘r’ which beyond this distance, the particle will not get trapped. Basically r is set to zero and then add a small increment each time you check a new initial condition – if the particle is captured (d-particle < 10^-10 units away form trap), then you increment farther out. This will continue until this condition is not true, and the graph of this specific instance will be shown. The problem is, it is not doing what I want it to, and I think there is something wrong with my syntax, not logic – I am used to other programming languages so I went about it like those, but I realize Wolfram is different so I’m not sure where I’m going wrong. Thanks for the help and let me know if you have any questions!

4

I am concerned that there are some simple misunderstandings of syntax of Mathematica code, which may look similar to other programming languages but may be very very different.

I am also concerned that the very small and very large numbers you are using in your differential equation may result in nothing but roundoff errors from the calculation.

Please study this and look at the result of running the calculation.

r = 0;(*defines the critical radius variable*)
For[i = 149999/10^5, i <= 150100/10^5, i += 0000015/10^5,
 s = NDSolve[{Rationalize[x'[t] - (2.5*10^-25) ((1.2*10^31 ((9*10^-6) - (1.5*10^-3) t)(5*10^5 ((9*10^-6) -
   (1.5*10^-3) t) - 1000000 x[t] + 3)(5*10^5 ((9*10^-6) - (1.5*10^-3) t) + 1000000 x[t] - 3)/(1000000000000
   ((9*10^-6) - (1.5*10^-3) t)^2 + 1000000000000 x[t]^2 -6000000 x[t] + 9)^(7/2)) - ((3 ((9*10^-6) -
   (1.5*10^-3) t)^3-12((9*10^-6)-(1.5*10^-3)t) x[t]^2))/(x[t]^2+((9*10^-6)-(1.5*10^-3) t)^2)^(7/2)),10^-25] == 0,
   x[0] == i*10^-6}, x, {t, 0, 1}][[1]];(*defines the differential equation to be solved with init.conditions*)
 Print["x[.006]=", x[t] /. s /. t -> 006/10^3];
 Print[Plot[Evaluate[x[t] /. s], {t, 0, 1/10}]];
 If[(x[t] /. s /. t -> 006/100) <= 10^-10,
  Print["True"]; r = r + 0000015/10^6,
  (*else*)
  Print["False"]; Print[Plot[Evaluate[x[t] /. s], {t, 0, 1/10}]],
  (*else*)
  Print["Neither True Nor False"]
  ];(*checks if zero at specified time-adds .000015 to r if true*)
 Print[r]
 ]

This only does a small number of iterations and prints some diagnostic information for each iteration to try to help you see what is being done.

I have made a number of changes to your code and I hope that these are correct. There was a ] in the wrong place in your For. This is one example of how Mathematica syntax and semantics are very different from other programming languages.

I cannot tell if your final If was intended to both increment r and display a plot when the calculation was <= 10^ -10 or if it meant to increment when true and Plot when false. If both were intended when true then you should put a semicolon between the increment and the Plot. The differing meaning of semicolon and comma is another example of how Mathematica differs.

Because of the very large and small numbers, I believe the result returned from NDSolve is often incorrect. You can look at the output plot and then try NDSolve[…, x, {t,0,1}, WorkingPrecision->50] instead of just NDSolve[…, x, {t,0,1}]. That will attempt to do all calculations to 50 significant digits instead of 16 digits. You should see that the results are now completely different. If you then change the precision to 100 you should see that the results are again completely different. But even with these there are a stream of warning and error messages from NDSolve. You may see that I introduced Rationalize[…, 10^-25] around your differential equation and made other changes to try to force all coefficients to be much more precise and avoid many of your cofficients having only one or two digits of precision as entered. I still do not believe that the result is well posed and open to accurate calculations yet. The warnings about the InterpolatingFunction values lying outside the range are another hint that there are still more things that need to be changed before these results can be used.

If you can compare character by character the code you posted and the revision I posted and try to understand what the changes were then perhaps you can discover any errors that I made and perhaps begin to track down what further changes need to be made.