Homework Help: If someone gives me the rundown on the 4th-order runge-kutta method, I'll love you

1. Apr 29, 2005

schattenjaeger

forever!

I missed a day of notes, I know for the second order

Y(n+1) =(approx.) Y(n)+K2, and I have the algorithm for finding k1 and then k2, how does this differ from the 4th order?

2. Apr 29, 2005

Hippo

A Fourth-Order Runge-Kutta Method

$$y_{n+1} = y_{n} + \frac{1}{6} \left( k_{1} + 2k_{2} + 2k_{3} + k_{4} \right)$$

where

$$k_{1} = h f \left( x_{n}, y_{n} \right)$$

$$k_{2} = h f \left( x_{n} + \frac{1}{2}h, y_{n} + \frac{1}{2}k_{1} \right)$$

$$k_{3} = h f \left( x_{n} + \frac{1}{2}h, y_{n} + \frac{1}{2}k_{2} \right)$$

$$k_{4} = h f \left( x_{n} + h, y_{n} + k_{3} \right)$$

Will you love me... forever? :!!)

Last edited by a moderator: Apr 29, 2005
3. Apr 29, 2005

schattenjaeger

Thanks! The advantage of that is it's more accurate, I assume?

Edit: FOREVARRRRRR

Last edited: Apr 29, 2005
4. Apr 29, 2005

Hippo

Yes.

5. Apr 29, 2005

Theelectricchild

6. Apr 29, 2005

schattenjaeger

Oh goody, now I'm having trouble using it right...
http://pacific.uta.edu/~qiming/Project_3.htm [Broken]

so starting with dx/dt = v and dv/dt = -kx/m - 10v for the diff. eq, is that part right at least?

then to find the position, I need k1=h*f(T(n),(Yn)), right?

so will that be...0 since initial velocity is 0?

then I can only assume for the velocity one I need

k1= h * -kx/m - 10v, initial velocity is 0, so k1 is just -kx/m?

Lordy this is confusing EDIT again: wait wait, k1=-kx/m - 10v, right? Because you need v in the final step so...um...help?

Edit: I guess I don't even need to use Runge-Kutta for dx/dt, do I? I just do it for dV/dt then use that to find the velocity at a bunch of points, and from that I can find the distance? Or...something

Last edited by a moderator: May 2, 2017
7. Apr 29, 2005

schattenjaeger

dangit, I can't figure out the position part....

I knoooow someone here knows how to do this:(

8. Apr 29, 2005

juvenal

Your equation for dv/dt is wrong since you are equating acceleration and force.

I'm not sure what you're confused about. The RK method is "plug-and-chug". There's not really any thinking involved. I do see that Hippo's notation might be a bit confusing, though. Instead of y in his equation, use a vector(x), and instead of x in his equation you want t. The k's are vector(k)'s. The subscripts in the k's above do not refer to time, btw, but to the intermediate steps in the RK to get from t_n to t_(n+1).

In your case, what you want is vector(x) = (x,v). Don't worry about the t since there is no explicit time dependence in the two differential eqns.

Now, start with (x_0, v_0). Use the RK method to get (x_1, v_1). Then use (x_1,v_1) to get (x_2,v_2). So on.

If you are still confused, please write out your work, step-by-step, getting from (x_0,v_0) to (x_1,v_1) explicitly, with as little extra verbiage as possible.

Last edited: Apr 29, 2005
9. Apr 29, 2005

schattenjaeger

Ok, thanks! I'm using the second-order RK, btw...
for the actual differential equation...

A=F/m, A=d^2x/dt^2
introduce V=dx/dt

dx/dt=V

dv/dt=(-kx-10v)/m -did I fix that right?
Xo=.08, Vo = 0
so I have
dx/dt=v
dv/dt=(-kx-10v)/m

To solve for x(t)
k1=h*(v)
k2=h*.08+(h*v/2)

and here's where I see my problem, x should be the variable, not v, 'cuz I'm going to need to find X(n+1) and plug that into the Xn+k2 equation to find step after step, right?

I'll have the same problem when I try to solve V(t), won't I?

Last edited: Apr 29, 2005
10. Apr 29, 2005

juvenal

I would write it out in vector notation (x,v). It's much easier to see that way.

As you said: vector(f(x,v)) = (v, (-kx - 10v)/m).

vector(k1) = h* (v_0, (-k*x_0 - 10*v_0)/m)

Now you need vector(k2). In order to calculate vector(k2), you need:

vector(x_0) + vector(k1)/2 =

(x_0,v_0) + (0.5 *h) * (v_0, (-k*x_0 - 10*v_0)/m) = Something simplified, call it "vector(a)" = (a1,a2)

Now, you need to calculate vector(f(vector(a))). That's equal to:

(a2, (-ka1 - 10a2)/m)

You know what a1, a2 are in terms of x0,v0, so you can rewrite f(vector(a)) in terms of x0, v0.

Then vector(k2) = h * f(vector(a)).

I hope you can take it from there.

Last edited: Apr 29, 2005
11. Apr 29, 2005

schattenjaeger

Thanks! I was just having trouble wrapping my brain around all that for some reason. The vector notation also makes it simpler to program, thanks!

12. Apr 29, 2005

The big advantage is that it is very accurate for an explicit method. It only requires function evaluations, no jacobians, so it's simple to implement.

13. Apr 29, 2005

schattenjaeger

Well then, the programming bit is annoying me now.

I'm running a pretty little do loop, and it looks like velocity is coming out right, but I'm starting the position at .08 and decreasing it by K2(since it's approaching x=0, or should be)and because k2 depends on x as well, the whole thing basically will never reach 0. That's a helluva spring, and a bit of an impossible one, sooo, the best thing I can figure is to do .08+k2, and treat .16 as the equilibrium. I think I may've had it right, I had the spring, starting at .08, with that friction force of 10v accelerate to the equilibrium, where it almost stopped cold and then just oscillates a little bit back and forth(a very little bit!)and the velocity does the same(accelerates to equilibrium, then oscillates)

does that sound reasonable?