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

Mathematica problem: two nigh-identical codes

  1. Mar 27, 2012 #1
    Hello, everyone.

    I have a rather bizzarre problem. I wrote a code that takes a matrix U[1] and calculates a matrix U[2] from it using simple math. I wrote three For loops, one to go across the U[n], the other two to read the matrix U[n] across all the lines and columns. I do it twice because it's a predictor/corrector kind of computaion. Pretty straightforward. Here's the code:

    For[n = 1, n <= 10, n++, UU2 = ConstantArray[0, {dim, dim}];
    UU3 = ConstantArray[0, {dim, dim}];
    UU = U[n];
    UU2[[1, All]] = UU[[1, All]]; UU2[[10, All]] = UU[[10, All]];
    UU2[[All, 1]] = UU[[All, 1]]; UU2[[All, 10]] = UU[[All, 10]];
    UU3 = UU2;
    For[i = 2, i <= dim - 1, i++,
    For[j = 2, j <= dim - 1, j++,
    UU2[[i, j]] =
    UU[[i, j]] - dt/dx*(G[UU[[i, j]]] - G[UU[[i - 1, j]]]) -
    dt/dx*(H[UU[[i, j]]] - H[UU[[i, j - 1]]]) + dt*S[UU[[i, j]]];
    ]
    ]
    For[i = 2, i <= dim - 1, i++,
    For[j = 2, j <= dim - 1, j++,
    UU3[[i, j]] =
    0.5*(UU[[i, j]] + UU2[[i, j]] -
    dt/dx*(G[UU2[[i + 1, j]]] - G[UU2[[i, j]]]) -
    dt/dx*(H[UU2[[i, j + 1]]] - H[UU2[[i, j]]]) +
    dt*S[UU2[[i, j]]]);
    ]
    ]
    U[n + 1] = UU3;
    ]


    The code above works. It gives me exactly what I need but it's suppose to be a bit unstable so I'm trying to write a different code now. The following code is nigh identical but for some reason Mathematica gives me all kinds of errors and fails to compute. Here's the faulty code:

    For[n = 1, n <= 3, n++,
    UU2 = ConstantArray[0, {dim, dim}];
    UU3 = ConstantArray[0, {dim, dim}];
    UU = U[n];
    UU2[[1, All]] = UU[[1, All]]; UU2[[dim, All]] = UU[[dim, All]];
    UU2[[All, 1]] = UU[[All, 1]]; UU2[[All, dim]] = UU[[All, dim]];
    UU3 = UU2;
    For[i = 2, i <= dim - 1, i++,
    For[j = 2, j <= dim - 1, j++,
    UU2[[i, j]] = UU[[i, j]] - dt/dx/2*(G[UU[[i, j]]] - G[UU[[i - 1, j]]]) + dt/2*S[UU[[i, j]]];
    ]
    ]
    For[i = 2, i <= dim - 1, i++,
    For[j = 2, j <= dim - 1, j++,
    UU3[[i, j]] = 0.5*(UU[[i, j]] - UU2[[i, j]] - dt/2/dx*(G[UU2[[i + 1, j]]] - G[UU2[[i, j]]]) +
    dt/2*S[UU2[[i, j]]]);
    ]
    ]
    U[n + 1] = UU3;
    ]

    The only difference I can see is that the calculation done in each step of the For loop is a bit shorter since it doesn't employ the H[UU2[[i,j]]] parts. I don't see how this could affect the overall scheme.

    These are the errors Mathematica reports when I execute the second piece of code:

    Set::write: Tag Times in Null Null U[2] is Protected. >>
    Part::partw: Part All of U[2] does not exist. >>
    Part::partd: Part specification U[2][[All,1]] is longer than depth of object. >>
    Part::partd: Part specification U[2][[All,10]] is longer than depth of object. >>
    Part::partw: Part 2 of U[2] does not exist. >>
    Part::partw: Part 2 of U[2] does not exist. >>
    General::stop: Further output of Part::partw will be suppressed during this calculation. >>
    Part::partd: Part specification U[2][[1,2]] is longer than depth of object. >>
    General::stop: Further output of Part::partd will be suppressed during this calculation. >>
    Set::write: Tag Times in Null Null U[3] is Protected. >>
    Set::write: Tag Times in Null Null U[4] is Protected. >>
    General::stop: Further output of Set::write will be suppressed during this calculation. >>

    Of course I clear all variables before running the second code and reload the definitions for G, H and S which are being used in the computation.
     
  2. jcsd
  3. Mar 27, 2012 #2
    This
    Set::write: Tag Times in Null Null U[2] is Protected.
    probably means at least in part you are missing a couple of semicolons. Yes, sometimes Mathematica will let you get away with that, and sometimes it will bite you. Sometimes without a semicolon Mathematica thinks two things adjacent to each other, like i j, means you want to multiply i and j. But it can also think that For[...] x means you want to multiply the Null returned from For by x.

    Stick semicolons between these two lines

    For[i=2,i<=dim,i++...]]
    For[i = 2, i <= dim - 1, i++,

    and between these two lines

    ]
    U[n + 1] = UU3;

    that is at the bottom of your code.

    That isn't going to fix all your problems and in particular not fix your "Part 2 of U[2] does not exist" problem, but it should get rid of all your "Tag Times in" problem.

    Now you haven't shown the rest of your code where U[] is defined. Perhaps you mean that U is a matrix and not a function. If my guessing is correct then U[[]] instead of U[] may make some of your problems go away. Maybe U[] really is a function, but then U[n+1]=UU3 looks suspicious. I can't tell what you really meant.
     
  4. Mar 28, 2012 #3
    SEMICOLONS!!!

    Blasted semicolons. That's all it took, at the ending of every ]. Thanks for the help, Bill. Sometimes I am such a ***.
     
  5. Mar 30, 2012 #4
    New problem, same code. Best not make another thread, I think.

    I used the codes mentioned before in a predictor-corrector scheme. This calls the same method four times, each time feeding it the new result.

    This is what the single method looks like:

    For[i = 2, i <= dim - 1, i++,
    For[j = 2, j <= dim - 1, j++,
    (*Izračun prediktorja Lx1*)
    UU2[[i, j]] =
    UU[[i, j]] -
    Dxy[UU[[i, j]]]*dt/dx/2*(G[UU[[i, j]]] - G[UU[[i - 1, j]]]) +
    Dxy[UU[[i, j]]]*dt/2*S[UU[[i, j]]];
    ];
    ];
    For[i = 2, i <= dim - 1, i++,
    For[j = 2, j <= dim - 1, j++,
    (*Izračun korektorja Lx1*)
    UU3[[i, j]] =
    0.5*(UU[[i, j]] - UU2[[i, j]] -
    Dxy[UU[[i, j]]]*
    dt/2/dx*(G[UU2[[i + 1, j]]] - G[UU2[[i, j]]]) +
    Dxy[UU[[i, j]]]*dt/2*S[UU2[[i, j]]]);
    ];
    ];

    I want to make that into a method or module (not sure how Mathematica calls it) so I can do something like this:
    UU3 = method[UU, UU2];

    (I used to use Matlab way more often than Mathematica so I'm still thinking in M-file ways.)

    I tried using the := way of identifying my methods, but it only seems to work for one input parameter. The corrector part uses both UU and UU2 to get UU3. Somehow writing it method[UU, UU2]:=... doesn't work.

    predCorr[UU_, UU2_] :=
    UU2 = ConstantArray[0, {dim, dim}];
    UU3 = ConstantArray[0, {dim, dim}];
    UU2[[1, All]] = UU[[1, All]]; UU2[[dim, All]] = UU[[dim, All]];
    UU2[[All, 1]] = UU[[All, 1]]; UU2[[All, dim]] = UU[[All, dim]];
    UU3 = UU2;
    For[i = 2, i <= dim - 1, i++,
    For[j = 2, j <= dim - 1, j++,
    (*Izračun prediktorja Lx1*)
    UU2[[i, j]] =
    UU[[i, j]] -
    Dxy[UU[[i, j]]]*dt/dx/2*(G[UU[[i, j]]] - G[UU[[i - 1, j]]]) +
    Dxy[UU[[i, j]]]*dt/2*S[UU[[i, j]]];
    ];
    ];
    For[i = 2, i <= dim - 1, i++,
    For[j = 2, j <= dim - 1, j++,
    (*Izračun korektorja Lx1*)
    UU3[[i, j]] =
    0.5*(UU[[i, j]] - UU2[[i, j]] -
    Dxy[UU[[i, j]]]*
    dt/2/dx*(G[UU2[[i + 1, j]]] - G[UU2[[i, j]]]) +
    Dxy[UU[[i, j]]]*dt/2*S[UU2[[i, j]]]);
    ];
    ];

    I get all sorts of errors, particularly these:

    Part::partd: Part specification UU[[1,All]] is longer than depth of object. >>
    Set::noval: Symbol UU2 in part assignment does not have an immediate value. >>
    Part::partd: Part specification UU[[10,All]] is longer than depth of object. >>
    Set::noval: Symbol UU2 in part assignment does not have an immediate value. >>
    Symbol::argx: Symbol called with 0 arguments; 1 argument is expected. >>
    Set::noval: Symbol UU2 in part assignment does not have an immediate value. >>
    Symbol::argx: Symbol called with 0 arguments; 1 argument is expected. >>
    Set::noval: Symbol UU2 in part assignment does not have an immediate value. >>
    Part::partd: Part specification UU[[2,2]] is longer than depth of object. >>
    Part::partd: Part specification UU[[2,2]] is longer than depth of object. >>

    I can tell that I can't simply copy/paste my working code into a delayed value scheme but I don't have a clue how to make it work.

    I made sure all semicolons are there this time. ;)
     
  6. Mar 30, 2012 #5
    Because of precedence I believe you are getting

    predCorr[UU_, UU2_] := UU2 = ConstantArray[0, {dim, dim}];
    (*and that is the end of your function definition and then *)
    UU3 = ConstantArray[0, {dim, dim}];
    (*and that is unrelated to your attempt to create a module *)

    It is difficult for me to guess what you want given your description.
    Perhaps this is something like what you are trying to do.

    mymodule[UU_] := Module[{UU2, UU3},
    For[i=2, i<=dim - 1, i++,
    For[j=2, j<=dim-1, j++,(*Izračun prediktorja Lx1*)
    UU2[[i,j]] = UU[[i,j]] - Dxy[UU[[i,j]]]* dt/dx/2* (G[UU[[i,j]]] - G[UU[[i-1,j]]]) + Dxy[UU[[i,j]]]*dt/2*S[UU[[i,j]]];];];
    For[i=2, i<=dim-1, i++,
    For[j=2, j<=dim-1, j++,(*Izračun korektorja Lx1*)
    UU3[[i,j]]=0.5*(UU[[i,j]] - UU2[[i,j]] - Dxy[UU[[i,j]]]*dt/2/dx*(G[UU2[[i+1,j]]] - G[UU2[[i,j]]]) + Dxy[UU[[i,j]]]*dt/2*S[UU2[[i,j]]]);];];
    UU3];
    mymodule[someUU]
     
    Last edited: Mar 30, 2012
  7. Mar 31, 2012 #6
    I made a mistake and left out the initialization of UU2 and UU3

    mymodule[UU_] := Module[{UU2, UU3},
    UU2 = ConstantArray[0, {dim, dim}]; (*Insert this *)
    UU3 = ConstantArray[0, {dim, dim}]; (* and this *)
    For[i=2, i<=dim - 1, i++,
    ...

    Sorry
     
  8. Apr 1, 2012 #7
    No problem. I was just going to ask about that. ;)
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Mathematica problem: two nigh-identical codes
Loading...