Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Mathematica loop problem.

  1. Nov 27, 2008 #1
    Hi there!

    I am having problem setting up a loop to run my program.

    Firstly, i have a DE solved using NDSolve.

    solu = NDSolve[{y''[x] == Sin[y[x]], y'[0] == 0, y[0] == a},
    y, {x, 0, 2}]

    and then to obtain my desired result, i need the initial condition (y[0]==a) to fulfill

    (y'[1] /. solu) + ((y[1] /. solu)/3)=0.0001

    then I will be given my result, <a>.

    I have manually found the result by changing the <a>. But I can’t always depend on this manual method as it is time consuming.

    Now I wish to just implement a loop that will automatically generate for <solu> from let say 0 to 1 with 0.01 interval.

    I’ve tried using <Catch> <Throw> but it fails. Here’s my code

    Catch[
    Do[
    Print[
    solu = NDSolve[{y''[x] == Sin[y[x]], y'[0] == 0, y[0] == a},
    y, {x, 0, 2}]];
    If[(y'[1] /. solu) + ((y[1] /. solu)/3) < 0.0001, Throw[a]], {a, 0, 1,
    0.1}]]


    But it fails. It just show a list of 10 NDSolve solutions without the <a> thrown out.
    I try again with <For> but it is even worse, I can’t even let the thing run.
    Please help me with this. Thank you very much.
     
  2. jcsd
  3. Nov 28, 2008 #2

    alphysicist

    User Avatar
    Homework Helper

    Hi DarkForXe,

    I believe this has an invalid comparison in your IF statement. You are checking whether:

    (y'[1] /. solu) + ((y[1] /. solu)/3) < 0.0001

    but mathematica does not know what to do with this. To see why, evaluate this on a line all by itself:

    (y'[1] /. solu) + ((y[1] /. solu)/3)

    and you'll get something like:

    {0.940501}

    and the curly brackets mean it is a list.

    What you need to do is get the value from the list, and then compare that value with 0.0001; there are several ways to do it. One way is to use the First[] function, which returns the first element of the list (makes sense), which works fine since your list only has one element.

    So if your code is changed to:

    Code (Text):

    Catch[
    Do[
    Print[
       solu = NDSolve[{y''[x] == Sin[y[x]], y'[0] == 0, y[0] == a},
      y, {x, 0, 2}]];
      If[ First[(y'[1] /. solu) + ((y[1] /. solu)/3)] < 0.0001, Throw[a]], {a, 0, 1,
       0.1}]]
     
    it throws the a=0 value for me.


    (and like I said, there are other ways to get it besides the First function.)
     
  4. Nov 28, 2008 #3
    Thanks alphysicist!!

    I do know that curly brackets means a list or and array of data but i do not know that it actually means number with very long decimal places and needed to be call out using the First function. I just assume that it is just a "normal" data and i didnt care much about the curly brackets.
    Besides this, i didn't know about the existance of First function as well.
    Guess i still have lots more to learn.

    Thanks again.:smile:
     
  5. Nov 28, 2008 #4

    alphysicist

    User Avatar
    Homework Helper

    Glad to help!

    And I don't believe it has anything to do with the number of decimal places in the number; it's just that mathematica does not compare a list with a number. If you give it this:

    0.9 < 0.1

    it knows what to do (return "false"), but with

    {0.9} < 0.1

    since the left and right side are different objects it will not compare them.




    Also, just off the top of my head, instead of the First function, you could have used the Last function, Part function, or just specified the element to be used, like this:

    If[ Last[ (y'[1] /. solu) + ((y[1] /. solu)/3) ] < 0.0001, Throw[a]]

    If[ Part[ (y'[1] /. solu) + ((y[1] /. solu)/3) ,1 ] < 0.0001, Throw[a]]

    If[ ( (y'[1] /. solu) + ((y[1] /. solu)/3) )[[1]] < 0.0001, Throw[a]]
     
  6. Nov 28, 2008 #5
    Thanks for clarifying. I am now at the final stage where i cleaning up the codes to make it as simple as possible now.

    Now i am having another simple problem yet i am clueless about it.
    Let say i have constant <a> which i declare earlier.
    Then i wish to set condition If a>0, it will go to process1, if a<0, it will go to process 2.

    I used the <If> command but it will just return me the answer after running both process. This happens also because i don't know how to assign a long code to a name, let say <pro1>?

    type pro1 and it will run the process.

    how to do it?
     
    Last edited: Nov 28, 2008
  7. Nov 28, 2008 #6

    alphysicist

    User Avatar
    Homework Helper

    I'm not sure what you are saying is happening here. Can you show some of the program that repeats the behavior (or if it's very long, could you write a toy program that shows the behavior)?
     
  8. Nov 30, 2008 #7
    Sorry, i have been away for the weekend and please pardon my broken english.
    Back to the question.
    Let say i have a <q>, let it be -3 in this case.
    so i in put

    q=-3

    then i wish to have a command that directs the program to go into s1 ONLY when <q> is less than 0. If It is more than 0, it will ONLY run s2.

    If i used
    If[q<0,s1,s2]
    It will just give me <s1> and run both process if i structure it as below.
    It is just a very simple thing but i really don't know how to build it out. Did i used the wrong thing and should be something else other than <If>?


    If[q<0,s1,s2];

    s1 = Catch[
    Do[solu =
    NDSolve[{y''[x] == -y[x] + y[x]^3, y'[0] == 0, y[0] == (1 + a)},
    y, {x, -2, 2}];
    If[First[(y'[1] /. solu) + ((y[1] /. solu)/(q))] > 0.0000001,
    Throw[a]], {a, 0, 1, 0.0001}]];
    s1=Plot[{Evaluate[y[x] /. solu]}, {x, -2, 2}]


    s2 = Catch[
    Do[solu =
    NDSolve[{y''[x] == -y[x] + y[x]^3, y'[0] == 0, y[0] == (1 - a)},
    y, {x, -2, 2}];
    If[First[(y'[1] /. solu) + ((y[1] /. solu)/(q))] < 0.0000001,
    Throw[a]], {a, 0, 1, 0.0001}]];
    Plot[{Evaluate[y[x] /. solu]}, {x, -2, 2}]
     
  9. Nov 30, 2008 #8

    alphysicist

    User Avatar
    Homework Helper


    Oh, I think I understand now! Here's the small changes I would do; let me know if it is what you are thinking about.

    First make s1 and s2 functions that you can call at a later time. So in a cell by themselves, put these (the red parts being the additions to your code):

    Code (Text):

    s1[COLOR="Red"][q_]:[/COLOR]= [COLOR="red"]{[/COLOR] Catch[
       Do[solu =
         NDSolve[{y''[x] == -y[x] + y[x]^3, y'[0] == 0, y[0] == (1 + a)},
          y, {x, -2, 2}];
        If[First[(y'[1] /. solu) + ((y[1] /. solu)/(q))] > 0.0000001,
         Throw[a]], {a, 0, 1, 0.0001}]];
    Plot[{Evaluate[y[x] /. solu]}, {x, -2, 2}] [COLOR="red"]}[/COLOR]


    s2[COLOR="red"][q_]:[/COLOR]= [COLOR="red"]{[/COLOR] Catch[
       Do[solu =
         NDSolve[{y''[x] == -y[x] + y[x]^3, y'[0] == 0, y[0] == (1 - a)},
          y, {x, -2, 2}];
        If[First[(y'[1] /. solu) + ((y[1] /. solu)/(q))] < 0.0000001,
         Throw[a]], {a, 0, 1, 0.0001}]];
    Plot[{Evaluate[y[x] /. solu]}, {x, -2, 2}] [COLOR="red"]}[/COLOR]
     
    Also notice that I removed the "s1=" from the beginning of line seven of s1; it looks like to me like that should not be there.

    (I guess you cannot just copy and paste these in mathematica or the red letters will get carried over and mess things up; but it's not much to add in by hand!)

    Then, later (after the functions have already been defined), you can call these functions by using, for example:

    Code (Text):


    q = 3;

    If [q<0 , s1[q], s2[q] ];

     

    Is that what you were looking for? (I wasn't sure if the plot was the only output you were wanting.)






    (By the way, your English seems fine to me!)
     
  10. Nov 30, 2008 #9
    It did solve part of my problem!

    I always stuck with this kind functions problem where i don't know how to input them in the memory and call them out.
    My real case is much more complicated and i will try my best to solve the problem before seeking for help again.

    Thanks again!
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: Mathematica loop problem.
  1. Mathematica looping (Replies: 3)

Loading...