How do I apply rk4 to a second order pde?

[DivergentSpectrum]

Im writing a program that calculates the trajectory of a particle in an arbitrary force field.
the force field is a vector function of position (x, y, z) AND velocity (x', y', z').
Rk4= runge kutta forth order method
Please help. Thanks!

 

[AlephZero]

Convert it into two coupled first-order equations.

If you have something like $m\frac{d^2x}{dt^2} = f(t,x,\frac{dx}{dt})$, take the two variables as $x$ and $v = \frac{dx}{dt}$.

The original equation then becomes the first-order equation $m\frac{dv}{dt} = f(t,x,v)$ and the other first-order equation is just $\frac{dx}{dt} = v$.

 

[DivergentSpectrum]

im not sure i understand... I am using a numerical method so once I've solved for velocity i don't have an equation to plug values into anymore, just a bunch of points. I am sure ill get it eventually though

 

[AlephZero]

If you have the RK algorithm that solves $\frac{dx}{dt} = f(x,t)$, it works exactly the same way when $x$ and $\frac{dx}{dt}$ are vectors of length 6, with components
$$x = \begin{bmatrix}u_x \\ u_y \\ u_z \\ v_x \\ v_y \\ v_z\end{bmatrix},\quad
\frac{dx}{dt} = \begin{bmatrix}\tfrac{du_x}{dt} \\ \tfrac{du_y}{dt} \\ \tfrac{du_z}{dt} \\
\tfrac{dv_x}{dt} \\ \tfrac{dv_y}{dt} \\ \tfrac{dv_z}{dt} \end{bmatrix}$$
Where $u_x$, $u_y$, $u_z$ are the three components of the position, and $v_x$, $v_y$, $v_z$ are the components of the velocity.

 

[DivergentSpectrum]

Thanks for the help. I think i have it figured out, unfortunately my 3d isn't very good so I am not sure if I am having a problem with graphics or if i did the algorithm wrong.

        while (i < n - 1)
        {
            k1vx = dt*ax(x[i], y[i], z[i], vx[i], vy[i], vz[i]);
            k1vy = dt*ax(x[i], y[i], z[i], vx[i], vy[i], vz[i]);
            k1vz = dt*ax(x[i], y[i], z[i], vx[i], vy[i], vz[i]);
    

            k1x = k1vx*dt;
            k1y=k1vy*dt;
            k1z=k1vz*dt;
   
            k2vx = dt*ax(x[i]+k1x/2, y[i]+k1y/2, z[i]+k1z/2, vx[i]+k1vx/2, vy[i]+k1vy/2, vz[i]+k1vz/2);
            k2vy = dt*ay(x[i]+k1x/2, y[i]+k1y/2, z[i]+k1z/2, vx[i] +k1vx / 2, vy[i] + k1vy / 2, vz[i] + k1vz / 2);
            k2vz = dt*az(x[i] + k1x / 2, y[i] + k1y / 2, z[i] + k1z / 2, vx[i] + dt * k1vx / 2, vy[i] + dt * k1vy / 2, vz[i] + dt * k1vz / 2);
       
            k2x=k1vx*dt;
            k2y=k1vy*dt;
            k2z = k1vz * dt;
         
            k3vx = dt*ax(x[i] + k2x / 2, y[i] + k2y / 2, z[i] + k2z / 2, vx[i] + k2vx / 2, vy[i] + k2vy / 2, vz[i] + k2vz / 2);
            k3vy = dt*ay(x[i] + k2x / 2, y[i] + k2y / 2, z[i] + k2z / 2, vx[i] + k2vx / 2, vy[i] + k2vy / 2, vz[i] + k2vz / 2);
            k3vz = dt*az(x[i] + k2x / 2, y[i] + k2y / 2, z[i] + k2z / 2, vx[i] + k2vx / 2, vy[i] + k2vy / 2, vz[i] + k2vz / 2);
         
            k3x = k3vx * dt;
            k3y = k3vy * dt;
            k3z = k3vz * dt;
       
            k4vx = dt*ax(x[i] + k3x, y[i] + k3y, z[i] + k3z, vx[i] + k3vx, vy[i] + k3vy, vz[i] + k3vz);
            k4vy = dt*ay(x[i] + k3x, y[i] + k3y, z[i] + k3z, vx[i] + k3vx, vy[i] + k3vy, vz[i] + k3vz);
            k4vz = dt*az(x[i] + k3x, y[i] + k3y, z[i] + k3z, vx[i] + k3vx, vy[i] + k3vy, vz[i] + k3vz);
        
            k4x = k4vx * dt;
             k4y=k4vy*dt;
             k4z = k4vz * dt;
      
            vx[i + 1] = vx[i] + (k1vx / 6 + k2vx / 3 + k3vx / 3 + k4vx);
            vy[i + 1] = vy[i] + (k1vy / 6 + k2vy / 3 + k3vy / 3 + k4vy);
            vz[i + 1] = vz[i] + (k1vz / 6 + k2vz / 3 + k3vz / 3 + k4vz);
       
           x[i + 1] = x[i] + (k1x / 6 + k2x / 3 + k3x / 3 + k4x / 6);
           y[i + 1] = y[i] + (k1y / 6 + k2y / 3 + k3y / 3 + k4y/ 6);
           z[i + 1] = z[i] + (k1z / 6 + k2z / 3 + k3z / 3 + k4z / 6);
            
            i++;
        }

 

[DivergentSpectrum]

it still won't work. eulers method works however, and i have some fascinating pictures:

elipseperihelion.jpeg shows an elipse due to gravity(inverse square), the advance of the perihelion is due to innaccuracies of the method (this is Newtonian gravity)
perfectelipse.jpeg shows the same orbit, only with more steps used.

Magnetic.jpeg shows the trajectory of a particle in a gravitational(or electrostatic) field and constant magnetic field pointing up along the z axis

imaginaryforce.jpeg shows the trajectory of a particle due to a central force that is proportional to 1/R. i believe part of the advancement of the perihelion is due to the the actual solution. (there should be 3 "petals")

imaginaryforce2.jpeg shows the trajectory of a particle due to a constant central force(Force strength doesn't depend on distance, but is always radial.)

red is the x axis, green is the y axis, and blue is the z axis

they kinda remind me of these things i had when i was a kid that draws similar designs.the 3d library i use is called helix3d for wpf. its really pretty easy to work with once you get over the learning curve

edit: problem solved. i also added mass consideration to the program, and a work energy calculation(using simpsons rule). if anyone wants to see the full source I am posting it in computers and technology.

 

[lpetrich]

I suggest trying your solver on some simple problem where the solution is relatively simple, like the 1D harmonic-oscillator problem. It should be easy to set up its force law and plug it in. As it runs, you may get some clues as to what is going wrong.

 

[Redbelly98]

FYI, this is not a partial differential equation. There is only one dependent variable, time.

I suggest trying your solver on some simple problem where the solution is relatively simple, like the 1D harmonic-oscillator problem. It should be easy to set up its force law and plug it in. As it runs, you may get some clues as to what is going wrong.

I agree. A classic problem with an exact solution to compare the numerical results. And it's 1-d, so the OP gets practice using two 1st order equations to model a single second-order equation without the complication of extra dimensions.

 

[lpetrich]

it still won't work. eulers method works however, ...

I suspect GIGO, garbage in, garbage out, because I don't see how RK4 would fail while Euler's method succeeds. Euler's method is the simplest method possible, and thus the dumbest.

I suggest that you take a closer look at your coding of RK4 and Euler's method. You ought to do it in a generic sort of fashion, with a generic vector of data and a generic function to operate on it. You'll find it much easier to try out on test cases that way. I don't know how you've been coding it, but in C++, here is what I would do. Create a generic RK4 class with a constructor that accepts the array length. The constructor will then allocate all the temporary space that it will use, and the destructor will then deallocate it. Or you can use some STL container class like "vector", a class which will do the deallocation for you. Also give it a virtual function that accepts as args the independent values, the data values, and the resulting function values, the last two as arrays.

To use this implementation, you create a subclass that gives both the function and the array-length value.

I'll now assess Euler's method and RK4 with a simple problem: y' = y with y(0) = 1 and going one step with size h. Its solution is
e^h ~ 1 + h + (1/2)*h² + (1/6)*h³ + (1/24)*h⁴ + (1/120)*h⁵ + ...

Euler's method:
f1 = 1
y = 1 + h

Midpoint method:
f1 = 1
f2 = 1 + (1/2)*h
y = 1 + h + (1/2)*h²

Heun's trapezoidal method:
f1 = 1
f2 = 1 + h
y = 1 + h + (1/2)*h²

RK4:
f1 = 1
f2 = 1 + (1/2)*h
f3 = 1 + (1/2)*h + (1/4)*h²
f4 = 1 + h + (1/2)*h² + (1/4)*h³
y = 1 + h + (1/2)*h² + (1/6)*h³ + (1/24)*h⁴

As you can see, these methods make the first few terms of the series for e^h. You can also see why RK4 is MUCH better than Euler's method.

 

[DivergentSpectrum]

actually i figured it out a while ago, thanks for all the help though. i saw this youtube video and apparently this is the proper method(for second order ode which i generalized to second order pde).

given x''=f(t,x,v) find x(t)
where v= x'

dx₁=h*v₀
dv₁=h*f(t₀,x₀,v₀)
dx₂=h*(v₀+dv₁/2)
dv₂=h*f(t₀+dt/2,x₀+dx₁/2,v₀+dv₁/2)
dx₃=h*(v₀+dv₂/2)
dv₃=h*f(t₀+dt/2,x₀+dx₂/2,v₀+dv₁/2)
dx₁=h*(v₀+dt)
dv₁=h*f(t₀+dt,x₀+dx3,v₀+dv1)

dx=(dx₁+2*dx₂+2*dx₃+dx₄)/6
dv=(dv₁+2*dv₂+2*dv₃+dv₄)/6
x₀+dx=x₁
v₀+dv=v₁

but i don't get why the function f is only evaluated using dv₁?

 

[lpetrich]

It looks like you made some mistakes. I will give the correct formulas for what I used: Euler, midpoint, trapezoid, and RK4.

We must solve dx/dt = f(t,x), with initial conditions t = t0 and x = x0 and advancing by h in t.

Euler's method:
f0 = f(t0, x0)
t = t0 + h
x = x0 + h*f0

Midpoint:
f0 = f(t0, x0)
f1 = f(t0 + h/2, x0 + (h/2)*f0)
t = t0 + h
x = x0 + h*f1

Trapezoid:
f0 = f(t0, x0)
f1 = f(t0 + h, x0 + h*f0)
t = t0 + h
x = x0 + (h/2)*(f0 + f1)

RK4:
f0 = f(t0, x0)
f1 = f(t0 + h/2, x0 + (h/2)*f0)
f2 = f(t0 + h/2, x0 + (h/2)*f1)
f3 = f(t0 + h, x0 + h*f2)
t = t0 + h
x = x0 + (h/6)*(f0 + 2*f1 + 2*f2 + f3)

In the general case, make x a vector and f return a vector of corresponding values.bluntwcrackrap, I have some questions about how you are writing your RK4 code.
What programming language are you using?
Are you using a generic array of variables x? Or are you letting them be specialized?
Are you using a generic function that accepts t and vector x and returns vector f(t,x)? Or are you writing some specialized function?
Have you tested your code on any problems with simple analytic solutions? Problems like the harmonic oscillator.

 

[DivergentSpectrum]

here is the video, don't get why it works, but yeah.
part of the problem searching for rk4 second order is that the search engine thinks you mean the runge kutta second order method lol figured it might help people also would be nice to know why it works this way.

Currently i have a problem adding a new feature to my program. (numerically calculating the force due to a current carrying wire given parametrically)
heres the thread:
https://www.physicsforums.com/threa...ectromagnetic-trajectory.771485/#post-4856175