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

Numerical second order pde solver

  1. Sep 1, 2014 #1
    Edit:whoops wrong forum mods please move

    2nd edit: I just had dinner then got back on the computer, input some points and saw a beautiful elipse.(complete with a fascinating flower petal design due to inaccuracies) Weird lol! No idea why it wasnt working before Now to implement RK4 bwahahaha

    im writing a program in WPF(an amalgam of c# and xaml) that tracks the position of a particle after a certain amount of time.
    the particle is placed in a force field, where force is a function of position and velocity.
    I used a crude numerical method, which seemed pretty obvious to me but for some reason it wont work. the particle appears to travel in a straight line as if there was no force.

    the number of increments and the final time to be calculated for is decided by the user, and the variable dt is calculated from that.
    vx,vy,vz are the velocity vector
    ax,ay,az are functions, and also the acceleration vector.

    heres the relevant code. im pretty sure the problem is in the first block.



    Code (Text):
          while (i < n - 1)
            {
       

                vx[i+1] = vx[i] + dt * ax(x[i], y[i], z[i], vx[i], vy[i], vz[i]);
                vy[i + 1] = vy[i] + dt * ay(x[i], y[i], z[i], vx[i], vy[i], vz[i]);
                vz[i + 1] = vz[i] + dt * az(x[i], y[i], z[i], vx[i], vy[i], vz[i]);


                x[i + 1] = x[i] + vx[i + 1] * dt;
                y[i + 1] = y[i] + vy[i + 1] * dt;
                z[i + 1] = z[i] + vz[i + 1] * dt;

               

                i++;
            }


    for verifiability, im using a gravitational field, also for simplicity i have eliminated force from the equation and im working with acceleration alone(ill add that later)
    here i have methods define the acceleration vectors:

    Code (Text):
         //acceleration vectors x componentvvvvvvv
            static double ax(double x, double y, double z, double vx, double vy, double vz)
            {
                double ax =  10000*(-x / Math.Pow(Math.Sqrt(x * x + y * y + z * z), 3));
                return ax;
            }
            //acceleration vectors x component^^^^^^


            //acceleration vectors y componentvvvvvvv
            static double ay(double x, double y, double z, double vx, double vy, double vz)
            {
                double ay =  10000*(-y / Math.Pow(Math.Sqrt(x * x + y * y + z * z), 3));
                return ay;
            }
            //acceleration vectors y component^^^^^^

            //acceleration vectors z componentvvvvvvv
            static double az(double x, double y, double z, double vx, double vy, double vz)
            {
                double az =  10000*(-z / Math.Pow(Math.Sqrt(x * x + y * y + z * z), 3));
                return az;

            }
       //acceleration vectors z component^^^^^^
     
     
    Last edited: Sep 1, 2014
  2. jcsd
  3. Sep 1, 2014 #2

    rcgldr

    User Avatar
    Homework Helper

    If you don't want to do RK4, you could try one trapezoidal corrector step:

    Code (Text):

        vx[i+1] = vx[i] + dt * ax(x[i], y[i], z[i], vx[i], vy[i], vz[i]);
        /* ... */
         x[i+1] =  x[i] + dt * vx[i];
        /* ... */
        vx[i+1] = vx[i] + .5 * dt * (ax(x[i], y[i], z[i], vx[i], vy[i], vz[i]) +
                  ax(x[i+1], y[i+1], z[i+1], vx[i+1], vy[i+1], vz[i+1]));
         x[i+1] =  x[i] + .5 *dt * (vx[i]  + vx[i+1]);
     
    update - see that you plan on using rk4, so no need to implement the corrector step method.
     
    Last edited: Sep 1, 2014
  4. Sep 1, 2014 #3

    AlephZero

    User Avatar
    Science Advisor
    Homework Helper

    The best way to find the problem this would be to print out all the variables for one time step (or inspect them with your IDE's debugger, if you can believe what it's telling you!) and check the values with a calculator.

    The code looks sensible, but if this is a snippet from a system that is built with tools that many people here have probably never used, the cause could be almost anything.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Numerical second order pde solver
  1. RK4 Solver Issue (Replies: 1)

  2. Sudoku Solver in C++ (Replies: 2)

Loading...