Understanding Second Order Algorithms

[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?

 

[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? :!)

 

[schattenjaeger]

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

Edit: FOREVARRRRRR

 

[Hippo]

Yes.

 

[Theelectricchild]

Whaaa! What about me!??

 

[schattenjaeger]

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

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

 

[schattenjaeger]

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

I knoooow someone here knows how to do this:(

 

[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.

 

[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?

 

[juvenal]

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 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.

 

[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!

 

[LeBrad]

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

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.

 

[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?